Libraries

# only use library paths in the anaconda environment

#.libPaths(grep('anaconda3', .libPaths(), value = T))

.libPaths()

## [1] "/data/users/cram/Project/Nyvac_096_Microbiome/packrat/lib/x86_64-pc-linux-gnu/3.6.0"
## [2] "/data/users/cram/Project/Nyvac_096_Microbiome/packrat/lib-ext"                      
## [3] "/data/users/cram/Project/Nyvac_096_Microbiome/packrat/lib-R"

R.version

##                _                                       
## platform       x86_64-pc-linux-gnu                     
## arch           x86_64                                  
## os             linux-gnu                               
## system         x86_64, linux-gnu                       
## status         beta                                    
## major          3                                       
## minor          6.0                                     
## year           2019                                    
## month          04                                      
## day            11                                      
## svn rev        76379                                   
## language       R                                       
## version.string R version 3.6.0 beta (2019-04-11 r76379)
## nickname       Planting of a Tree

sessionInfo()

## R version 3.6.0 beta (2019-04-11 r76379)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] compiler_3.6.0  backports_1.1.2 magrittr_1.5    rprojroot_1.3-2
##  [5] htmltools_0.3.6 tools_3.6.0     yaml_2.2.0      Rcpp_0.12.19   
##  [9] stringi_1.2.4   rmarkdown_1.10  knitr_1.20      stringr_1.3.1  
## [13] digest_0.6.18   packrat_0.4.8-1 evaluate_0.12

# https://stackoverflow.com/questions/46354826/have-a-function-that-calls-library-and-takes-either-a-package-or-its-name-as-inp


# Also return package version when loading in packages
# accept strings or functions
libver <- function(pac){

    pac <- as.character(substitute(pac))
    library(pac, character.only=TRUE)
    packageVersion(pac)
    }

#libver("dada2")
#libver("ggplot2")

libver("Cairo")

## [1] '1.5.9'

# Much of the data handling
libver('phyloseq')

## Registered S3 methods overwritten by 'ggplot2':
##   method         from 
##   [.quosures     rlang
##   c.quosures     rlang
##   print.quosures rlang

## Registered S3 method overwritten by 'dplyr':
##   method               from  
##   as.data.frame.tbl_df tibble

## [1] '1.22.3'

# A bunch of environments, including ggplot, dplyr, tidyr, and broom, which I use a lot
libver('tidyverse')

## Registered S3 method overwritten by 'rvest':
##   method            from
##   read_xml.response xml2

## ── Attaching packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──

## ✔ ggplot2 3.1.0     ✔ purrr   0.2.5
## ✔ tibble  1.4.2     ✔ dplyr   0.7.7
## ✔ tidyr   0.8.2     ✔ stringr 1.3.1
## ✔ readr   1.1.1     ✔ forcats 0.3.0

## ── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

## [1] '1.2.1'

# Mostly for concatenating ggplots
library(gridExtra); packageVersion("gridExtra")

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

## [1] '2.3'

# I use this surprisingly not a lot here.
library(vegan); packageVersion("vegan")

## Loading required package: permute

## Loading required package: lattice

## This is vegan 2.5-3

## [1] '2.5.3'

# For making trees
# libver('phangorn')

# A prerequesite to phangorn
# libver("DECIPHER")

# Some pre-processing stuff
# libver("dada2")

# I usually reshape with tidyverse tools now, but melt and cast are often easier in a pinch
# libver("reshape2")

# For replacing NaNs without too much thought.
# libver("imputeMissings")

# Deal with proportional data, especially useful for calculating proportionality phi
libver('compositions')

## Loading required package: tensorA

## 
## Attaching package: 'tensorA'

## The following object is masked from 'package:base':
## 
##     norm

## Loading required package: robustbase

## Loading required package: energy

## Loading required package: bayesm

## Welcome to compositions, a package for compositional data analysis.
## Find an intro with "? compositions"

## 
## Attaching package: 'compositions'

## The following objects are masked from 'package:stats':
## 
##     cor, cov, dist, var

## The following objects are masked from 'package:base':
## 
##     %*%, scale, scale.default

## [1] '1.40.2'

# Works with tidyverse to make model output tidy
libver('broom')

## [1] '0.5.0'

# Make pretty tables
libver('knitr')

## [1] '1.20'

libver('kableExtra')

## [1] '0.9.0'

# Let those pretty tables actually show up in a jupyter notebook
#library('IRdisplay')

# For bootstrapping
libver('boot')

## 
## Attaching package: 'boot'

## The following object is masked from 'package:robustbase':
## 
##     salinity

## The following object is masked from 'package:lattice':
## 
##     melanoma

## [1] '1.3.20'

# Calculate kernel regressions
libver("MiRKAT")

## Loading required package: survival

## 
## Attaching package: 'survival'

## The following object is masked from 'package:boot':
## 
##     aml

## The following object is masked from 'package:robustbase':
## 
##     heart

## Loading required package: PearsonDS

## Loading required package: GUniFrac

## Loading required package: ape

## 
## Attaching package: 'ape'

## The following object is masked from 'package:compositions':
## 
##     balance

## Loading required package: matrixStats

## 
## Attaching package: 'matrixStats'

## The following objects are masked from 'package:robustbase':
## 
##     colMedians, rowMedians

## The following object is masked from 'package:dplyr':
## 
##     count

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

## [1] '1.0.1'

libver("car")

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:boot':
## 
##     logit

## The following object is masked from 'package:dplyr':
## 
##     recode

## The following object is masked from 'package:purrr':
## 
##     some

## [1] '3.0.2'

#libver(mclust)

#libver(chemometrics)

libver(purrrlyr)

## [1] '0.0.3'

libver('qvalue')

## [1] '2.10.0'

libver("breakaway")

## [1] '4.6.11'

Functions

I have put the functions in a library file

source('libraries/library096.R')

Data

# Set upOriginal to false, if you want to used user-reprocessed data.
# Results may differ slightly from those in the manuscript due to inter-run variation
# especially in the tree-ing algorithm.
 upOriginal <- TRUE
# upOriginal <- FALSE

# For permutation tests, how fast do things need to run
# 9999 for most runs, 99999 for publication quality ones suggested
jnperm <- 99999

# Data paths
getwd()

## [1] "/data/users/cram/Project/Nyvac_096_Microbiome"

(mapping_file_path <- file.path('data', 'mapping_file_096a.csv'))

## [1] "data/mapping_file_096a.csv"

(immune_file_path <- file.path('data', 'immune096b.csv'))

## [1] "data/immune096b.csv"

if(upOriginal){
     seqtab_file_path <- file.path('data', 'seqtab.nochimNov2017.csv')
     taxa_file_path <- file.path('data', 'TaxaNov2017.csv')
     tree_path <- file.path('data', 'phylogeny096NovTree.tre')
    } else {
     seqtab_file_path <- file.path('data1', 'seqtab.nochimMar2018.csv')
     taxa_file_path <- file.path('data1', 'TaxaMar2018.csv')
     tree_path <- file.path('data1', 'phylogeny096Mar2018tre.tre')
}

seqtab_file_path

## [1] "data/seqtab.nochimNov2017.csv"

taxa_file_path

## [1] "data/TaxaNov2017.csv"

tree_path

## [1] "data/phylogeny096NovTree.tre"

# Sequence data
seqtab.nochim.data <- read.csv(seqtab_file_path)

seqtabNames = gsub('\\.', '-',
    gsub('.fastq', '', seqtab.nochim.data$X)
                   )

seqtab.nochim = as.matrix(seqtab.nochim.data[,-1])
rownames(seqtab.nochim) = seqtabNames

# Taxa names
taxa.data <- read.csv(taxa_file_path)
taxa = taxa.data[,-1]

## I reverse complemented the sequences to generate the taxonomy
# (but only in this latest re-run, not the original)
## The following undoes that reverse complement to get original sequence
#rownames(taxa) = dada2:::rc(taxa.data[,1]) 

if(upOriginal){
    rownames(taxa) = (taxa.data[,1])} else {
    rownames(taxa) = dada2:::rc(taxa.data[,1]) 
}

taxa <- as.matrix(taxa)

# Mapping file
mapping.data <- read_csv(mapping_file_path) %>%
mutate(pub_id = sapply(pub_id,  function(x) {as.numeric(gsub("096-", "", x))}))

## Parsed with column specification:
## cols(
##   SampleID = col_character(),
##   BarcodeSequence = col_character(),
##   LinkerPrimerSequence = col_character(),
##   ReversePrimer = col_character(),
##   run_prefix = col_character(),
##   pub_id = col_character(),
##   Sex = col_character(),
##   Visit = col_integer(),
##   visitRank = col_integer(),
##   RXCode = col_character(),
##   Description = col_character()
## )

## Warning: The `printer` argument is deprecated as of rlang 0.3.0.
## This warning is displayed once per session.

#mapping = mapping.data[,-1]
#rownames(mapping) = mapping.data[,1]
#mapping <- as.matrix(mapping)

head(mapping.data)

## # A tibble: 6 x 11
##   SampleID BarcodeSequence LinkerPrimerSeq… ReversePrimer run_prefix pub_id
##   <chr>    <chr>           <chr>            <chr>         <chr>       <dbl>
## 1 Sample-… TTACGC          GCGGACTACCVGGGT… GCACTCCTRCGG… JH3VWZ201     282
## 2 Sample-… TTAGGT          GCGGACTACCVGGGT… GCACTCCTRCGG… JH3VWZ201     282
## 3 Sample-… TTCCAC          GCGGACTACCVGGGT… GCACTCCTRCGG… JH3VWZ201     176
## 4 Sample-… TTGTAC          GCGGACTACCVGGGT… GCACTCCTRCGG… JH3VWZ201     176
## 5 Sample-… TTGTGT          GCGGACTACCVGGGT… GCACTCCTRCGG… JH3VWZ201     123
## 6 Sample-… TATCAC          GCGGACTACCVGGGT… GCACTCCTRCGG… JH3VWZ201      49
## # ... with 5 more variables: Sex <chr>, Visit <int>, visitRank <int>,
## #   RXCode <chr>, Description <chr>

# Immune Data
immune.data0 <- read_csv(immune_file_path)

## Parsed with column specification:
## cols(
##   visitno = col_integer(),
##   rx_code = col_character(),
##   type = col_character(),
##   antigen = col_character(),
##   mag = col_double(),
##   mag_bl = col_double(),
##   response = col_integer(),
##   day = col_integer(),
##   month = col_double(),
##   ct = col_character(),
##   response_j = col_double(),
##   assay = col_character(),
##   pub_id = col_character()
## )

immune.data <- mutate(immune.data0, pub_id = sapply(pub_id,  function(x) {as.numeric(gsub("096-", "", x))}))
levels(immune.data$antigen) <- gsub("[ /]", ".", levels(immune.data$antigen))

head(immune.data)

## # A tibble: 6 x 13
##   visitno rx_code type  antigen   mag mag_bl response   day month ct   
##     <int> <chr>   <chr> <chr>   <dbl>  <dbl>    <int> <int> <dbl> <chr>
## 1       5 CTRL    IgA   gp41     180.   214.        0    42   1.5 CTRL 
## 2       7 CTRL    IgA   gp41     175.   214.        0    98   3.5 CTRL 
## 3       9 CTRL    IgA   gp41     177.   214.        0   182   6.5 CTRL 
## 4      12 CTRL    IgA   gp41     173    214.        0   364  12   CTRL 
## 5       5 CTRL    IgA   p24      447.   454.        0    42   1.5 CTRL 
## 6       7 CTRL    IgA   p24      767.   454.        0    98   3.5 CTRL 
## # ... with 3 more variables: response_j <dbl>, assay <chr>, pub_id <dbl>

# Phylogenetic tree
seqs <- dada2::getSequences(seqtab.nochim)
names(seqs) <- seqs

pt <- ape::read.tree(file=tree_path)

pt2 <- phangorn::midpoint(pt)

immune.data$antigen %>% unique

## [1] "gp41"               "p24"                "Con.6.gp120.B"     
## [4] "ZM96.gp140"         "gp70_B.CaseA_V1_V2" "ANY.ENV.PTEG"

Save options to a variable

par0 <- options()

Initial Taxonomic filter.

Some investegation suggested by the phyloseq tutorials to identify phyla for removal, and to identify an abundance threshold

psN

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 960 taxa and 47 samples ]
## sample_data() Sample Data:       [ 47 samples by 6 sample variables ]
## tax_table()   Taxonomy Table:    [ 960 taxa by 8 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 960 tips and 959 internal nodes ]

psN %>% subset_samples(!is.na(pub_id))

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 960 taxa and 38 samples ]
## sample_data() Sample Data:       [ 38 samples by 6 sample variables ]
## tax_table()   Taxonomy Table:    [ 960 taxa by 8 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 960 tips and 959 internal nodes ]

# skip the blanks
psN %>% subset_samples(!is.na(pub_id)) %>%
# OTUs must be characterized to phylum
subset_taxa(!is.na(Phylum)& !Phylum %in% c("", "uncharacterized")) -> psN_hasPhylum
psN_hasPhylum

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 929 taxa and 38 samples ]
## sample_data() Sample Data:       [ 38 samples by 6 sample variables ]
## tax_table()   Taxonomy Table:    [ 929 taxa by 8 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 929 tips and 928 internal nodes ]

# from 960 to 929 otus

Identifying and removing phyla with very few taxa in them

prevdf = apply(X = otu_table(psN_hasPhylum),
                 MARGIN = ifelse(taxa_are_rows(psN_hasPhylum), yes = 1, no = 2),
                 FUN = function(x){sum(x > 0)})
# Add taxonomy and total read counts to this data.frame
prevdf = data.frame(Prevalence = prevdf,
                      TotalAbundance = taxa_sums(psN_hasPhylum),
                      tax_table(psN_hasPhylum))

plyr::ddply(prevdf, "Phylum", function(df1){cbind(mean(df1$Prevalence),sum(df1$Prevalence))})

##            Phylum        1    2
## 1  Actinobacteria 7.953846  517
## 2   Bacteroidetes 9.414634 1930
## 3   Elusimicrobia 2.000000    2
## 4      Firmicutes 9.711604 5691
## 5    Fusobacteria 5.875000   94
## 6  Proteobacteria 6.960000  348
## 7   Synergistetes 3.750000   15
## 8     Tenericutes 4.000000    4
## 9 Verrucomicrobia 6.000000    6

filterPhyla = c("Verrucomicrobia", "Tenericutes", "Elusimicrobia", "Synergistetes")
psN_MainPhyla = subset_taxa(psN_hasPhylum, !Phylum %in% filterPhyla)
psN_MainPhyla

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 922 taxa and 38 samples ]
## sample_data() Sample Data:       [ 38 samples by 6 sample variables ]
## tax_table()   Taxonomy Table:    [ 922 taxa by 8 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 922 tips and 921 internal nodes ]

# Determining abundance threshold
prevdf1 = subset(prevdf, Phylum %in% get_taxa_unique(psN_MainPhyla, "Phylum"))
ggplot(prevdf1, aes(TotalAbundance, Prevalence / nsamples(psN_hasPhylum),color=Phylum)) +
  # Include a guess for parameter
  geom_hline(yintercept = 0.05, alpha = 0.5, linetype = 2) + geom_point(size = 2, alpha = 0.7) +
  scale_x_log10() +  xlab("Total Abundance") + ylab("Prevalence [Frac. Samples]") +
  facet_wrap(~Phylum) + theme(legend.position="none")

## Warning: Transformation introduced infinite values in continuous x-axis

# Determining abundance threshold
prevdf1 = subset(prevdf, Phylum %in% get_taxa_unique(psN_MainPhyla, "Phylum"))
ggplot(prevdf1, aes(TotalAbundance, Prevalence / nsamples(psN_hasPhylum),color=Phylum)) +
  # Include a guess for parameter
  geom_hline(yintercept = 0.05, alpha = 0.5, linetype = 2) + geom_point(size = 2, alpha = 0.7) +
  scale_x_continuous() +  xlab("Total Abundance") + ylab("Prevalence [Frac. Samples]") +
  facet_wrap(~Phylum) + theme(legend.position="none")

Constructing psN2

(phyloseq object of relative abundances)

And psN1 (phyloseq object of counts)

psN1A and psN2A include all participants, and will be used to look at variability within participants

psN %>%
# add all the immune data
phylo_join(immune.table, by = "pub_id") %>%
# only use data from humans (no extraction controls)
subset_samples(is.finite(muVisitRank)) %>%
# only otus from known taxa that show up frequently enough
subset_taxa(!is.na(Phylum)& !Phylum %in% c("", "uncharacterized")) %>%
subset_taxa(!Phylum %in% filterPhyla) %>%
# only otus that show up in at least 10% of samples
prevalence_filter_taxa %>%
# convert to relative abundance

tag_phyloseq%>%
# Instead of naming each taxon with its full sequence, we use the "tag" instead
swap.phyloseq.taxnames %>%
pass -> psN1A # Save pre relative abundance transformation

## Warning: `list_len()` is deprecated as of rlang 0.2.0.
## Please use `new_list()` instead.
## This warning is displayed once per session.

## Warning: Setting row names on a tibble is deprecated.

# add is-male
manColumn <- psN1A %>% sample_data %>% as('data.frame') %>% rownames_to_column  %>% mutate(isMale = testIsMaleVec(sex)) %>% dplyr::select(rowname, isMale)
psN1A <- phylo_join(psN1A, manColumn, by = 'rowname')

## psN2 is like psN1 but with relative abundances
psN1A %>%
transform_sample_counts(function(x) {x/sum(x)}) %>%
# The "tag" is a new name that takes into account the rest of the taxonomy data
# the tag may need to be updated after any agglomeration
pass-> psN2A

# filter to just microbiome visit 1 and experemental treatments
psN1A %>%
subset_samples(muVisitRank == 1) %>%
subset_samples(ct == "T") %>%
pass -> psN1

psN2A %>%
subset_samples(muVisitRank == 1) %>%
subset_samples(ct == "T") %>%
pass -> psN2

# Calculate weighted unifrac distances and role those in.
psN2.wuf <- phyloseq::distance(psN2, method = "wunifrac")
psN2.pcoa <- capscale(psN2.wuf ~ 1)
psN2.pcoa.df <- psN2.pcoa %>% scores(display = "sites") %>%
        as.data.frame %>% 
        rownames_to_column %>% 
        dplyr::select('rowname', 'MDS1', 'MDS2') %>%
        mutate(rMDS1 = rank(MDS1)) %>% # rank order of MDS1
        mutate(rrMDS1 = formatC(format = "d", rMDS1, flag = "0", width=ceiling(log10(max(rMDS1))))) %>%
        unite(newname, rrMDS1, rowname, sep = "_", remove = FALSE) %>%
        dplyr::select(-rrMDS1)

psN2 %>%
phylo_join(
    psN2.pcoa.df,
    by = 'rowname'
) -> psN2

## Even if the data are counts, 
## the weighted unifrac pcoa is still done on the relative abundances
psN1 %>%
phylo_join(
    psN2.pcoa.df,
    by = 'rowname'
) -> psN1

psN2A

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 651 taxa and 38 samples ]
## sample_data() Sample Data:       [ 38 samples by 31 sample variables ]
## tax_table()   Taxonomy Table:    [ 651 taxa by 10 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 651 tips and 650 internal nodes ]

psN1A

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 651 taxa and 38 samples ]
## sample_data() Sample Data:       [ 38 samples by 31 sample variables ]
## tax_table()   Taxonomy Table:    [ 651 taxa by 10 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 651 tips and 650 internal nodes ]

psN2

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 651 taxa and 21 samples ]
## sample_data() Sample Data:       [ 21 samples by 35 sample variables ]
## tax_table()   Taxonomy Table:    [ 651 taxa by 10 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 651 tips and 650 internal nodes ]

psN1

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 651 taxa and 21 samples ]
## sample_data() Sample Data:       [ 21 samples by 35 sample variables ]
## tax_table()   Taxonomy Table:    [ 651 taxa by 10 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 651 tips and 650 internal nodes ]

Magnitude and variance of vaccine response per group

colnames(immune.data)

##  [1] "visitno"    "rx_code"    "type"       "antigen"    "mag"       
##  [6] "mag_bl"     "response"   "day"        "month"      "ct"        
## [11] "response_j" "assay"      "pub_id"

Plan. Calculate mean and variance of each antigen-type at visit 9, and at baseline.

#
geomean <- function(x, na.rm = FALSE, trim = 0, ...)
{
exp(mean(log(x, ...), na.rm = na.rm, trim = trim, ...))
}
 
geosd <- function(x, na.rm = FALSE, ...)
{
exp(sd(log(x, ...), na.rm = na.rm, ...))
}

immune.data %>% dplyr::filter(visitno %in% c(9), ct == "T") %>% dplyr::select(pub_id, antigen, visitno, mag, mag_bl, type) %>% group_by(antigen, type) %>%
summarize(mean_mag = geomean(mag), mean_bl = geomean(mag_bl), sd_mag = geosd(mag), sd_bl = geosd(mag_bl)) %>%
mutate(var_over_mean_mag = sd_mag^2/mean_mag, var_over_mean_bl = sd_bl^2/mean_bl)

## # A tibble: 8 x 8
## # Groups:   antigen [6]
##   antigen type  mean_mag mean_bl sd_mag sd_bl var_over_mean_m…
##   <chr>   <chr>    <dbl>   <dbl>  <dbl> <dbl>            <dbl>
## 1 ANY.EN… CD4+   9.20e-2   NA      2.51 NA           68.7     
## 2 Con.6.… IgG    1.14e+4    2.54   8.33  4.64         0.00607 
## 3 gp41    IgA    9.87e+1   15.9    6.42  7.04         0.417   
## 4 gp41    IgG    5.56e+3  208.     3.76  4.91         0.00254 
## 5 gp70_B… IgG    1.98e+3    1.48   7.40  2.66         0.0277  
## 6 p24     IgA    2.50e+3  422.     4.02  2.48         0.00646 
## 7 p24     IgG    2.55e+4  152.     2.44  6.34         0.000233
## 8 ZM96.g… IgG    1.81e+3    8.83   7.07  9.01         0.0276  
## # ... with 1 more variable: var_over_mean_bl <dbl>

Mean variance and variance over mean of each value.

immune.data %>% dplyr::filter(visitno %in% c(9), ct == "T") %>% dplyr::select(pub_id,rx_code, antigen, visitno, mag, mag_bl, type) %>% group_by(rx_code, antigen, type) %>%
summarize(mean_mag = geomean(mag), mean_bl = geomean(mag_bl), sd_mag = geosd(mag), sd_bl = geosd(mag_bl)) %>%
mutate(var_over_mean_mag = sd_mag^2/mean_mag, var_over_mean_bl = sd_bl^2/mean_bl) %>% 
gather(key = "meas", value = "value", mean_mag, mean_bl, sd_mag, sd_bl, var_over_mean_mag, var_over_mean_bl) %>% 
group_by(antigen, type, meas) %>% summarize(mean_val = mean(value)) %>%
spread(key = "meas", value = "mean_val")

## # A tibble: 8 x 8
## # Groups:   antigen, type [8]
##   antigen type  mean_bl mean_mag sd_bl sd_mag var_over_mean_bl
##   <chr>   <chr>   <dbl>    <dbl> <dbl>  <dbl>            <dbl>
## 1 ANY.EN… CD4+    NA     9.76e-2 NA      2.37          NA     
## 2 Con.6.… IgG      2.60  1.21e+4  4.77   8.59           9.04  
## 3 gp41    IgA     17.4   1.06e+2  7.01   6.62           4.30  
## 4 gp41    IgG    215.    6.43e+3  4.94   3.51           0.135 
## 5 gp70_B… IgG      1.51  1.99e+3  2.61   7.93           4.64  
## 6 p24     IgA    428.    3.34e+3  2.51   3.38           0.0147
## 7 p24     IgG    155.    2.56e+4  6.94   2.44           0.466 
## 8 ZM96.g… IgG      8.85  1.98e+3  9.47   7.15          10.6   
## # ... with 1 more variable: var_over_mean_mag <dbl>

As above, but this time, calculated seperately for each treatment group and then those values averaged. ZM96 and gp70 surprisingly large variance over mean. Maybe should look at delta magnitude

immune.data%>% head

## # A tibble: 6 x 13
##   visitno rx_code type  antigen   mag mag_bl response   day month ct   
##     <int> <chr>   <chr> <chr>   <dbl>  <dbl>    <int> <int> <dbl> <chr>
## 1       5 CTRL    IgA   gp41     180.   214.        0    42   1.5 CTRL 
## 2       7 CTRL    IgA   gp41     175.   214.        0    98   3.5 CTRL 
## 3       9 CTRL    IgA   gp41     177.   214.        0   182   6.5 CTRL 
## 4      12 CTRL    IgA   gp41     173    214.        0   364  12   CTRL 
## 5       5 CTRL    IgA   p24      447.   454.        0    42   1.5 CTRL 
## 6       7 CTRL    IgA   p24      767.   454.        0    98   3.5 CTRL 
## # ... with 3 more variables: response_j <dbl>, assay <chr>, pub_id <dbl>

immune.data %>% 
mutate(mag_delta = mag / mag_bl) %>%
mutate(mag_delta = if_else(mag_delta < 1, 1, mag_delta)) %>%
dplyr::filter(visitno %in% c(9), ct == "T") %>% dplyr::select(pub_id,rx_code, antigen, visitno, mag, mag_bl, mag_delta, type) %>% group_by(rx_code, antigen, type) %>%
summarize(mean_mag = geomean(mag), mean_bl = geomean(mag_bl), mean_delta = geomean(mag_delta), sd_mag = geosd(mag), sd_bl = geosd(mag_bl), sd_delta = geosd(mag_delta)) %>%
mutate(var_over_mean_mag = sd_mag^2/mean_mag, var_over_mean_bl = sd_bl^2/mean_bl, var_over_mean_delta = sd_delta^2/mean_delta) %>% 
gather(key = "meas", value = "value", mean_mag, mean_bl, sd_mag, sd_bl, var_over_mean_mag, var_over_mean_bl, mean_delta, sd_delta, var_over_mean_delta) %>% 
group_by(antigen, type, meas) %>% summarize(mean_val = mean(value)) %>%
spread(key = "meas", value = "mean_val") %>%
arrange(type)

## # A tibble: 8 x 11
## # Groups:   antigen, type [8]
##   antigen type  mean_bl mean_delta mean_mag sd_bl sd_delta sd_mag
##   <chr>   <chr>   <dbl>      <dbl>    <dbl> <dbl>    <dbl>  <dbl>
## 1 ANY.EN… CD4+    NA         NA     9.76e-2 NA       NA      2.37
## 2 gp41    IgA     17.4        8.97  1.06e+2  7.01     5.63   6.62
## 3 p24     IgA    428.         7.27  3.34e+3  2.51     3.89   3.38
## 4 Con.6.… IgG      2.60    5310.    1.21e+4  4.77    13.6    8.59
## 5 gp41    IgG    215.        35.8   6.43e+3  4.94     6.98   3.51
## 6 gp70_B… IgG      1.51    1417.    1.99e+3  2.61     9.60   7.93
## 7 p24     IgG    155.       176.    2.56e+4  6.94     7.87   2.44
## 8 ZM96.g… IgG      8.85     239.    1.98e+3  9.47    13.9    7.15
## # ... with 3 more variables: var_over_mean_bl <dbl>,
## #   var_over_mean_delta <dbl>, var_over_mean_mag <dbl>

Number of participants per group

immune.data %>% 
group_by(rx_code) %>%
summarize(Unique_ids = n_distinct(pub_id))

## # A tibble: 5 x 2
##   rx_code Unique_ids
##   <chr>        <int>
## 1 CTRL            16
## 2 T1              20
## 3 T2              20
## 4 T3              20
## 5 T4              20

donor.immune %>% 
group_by(rx_code) %>%
summarize(Unique_ids = n_distinct(pub_id))

## # A tibble: 5 x 2
##   rx_code Unique_ids
##   <chr>        <int>
## 1 CTRL             2
## 2 T1               7
## 3 T2               3
## 4 T3               6
## 5 T4               5

psN2 %>% sample_data %>% data.frame %>% group_by(muVisit) %>% summarize(n = length(muVisit))

## # A tibble: 3 x 2
##   muVisit     n
##     <int> <int>
## 1       2     3
## 2       9    11
## 3      12     7

psN2 %>% sample_data %>% data.frame %>% group_by(rx_code) %>% summarize(n = length(muVisit))

## # A tibble: 4 x 2
##   rx_code     n
##   <chr>   <int>
## 1 T1          7
## 2 T2          3
## 3 T3          6
## 4 T4          5

sampleBreakdown <- psN2 %>% sample_data %>% data.frame %>% group_by(muVisit, rx_code) %>% summarize(n = length(muVisit)) %>% spread(key = rx_code, value = n, fill = 0) %>% mutate(total = T1 + T2 + T3 + T4)

sbcs <- colSums(sampleBreakdown)

sampleBreakdown <- bind_rows(sampleBreakdown, sbcs)

sampleBreakdown %>%
kable("html", escape = F, digits = 3, align = 'c') %>%
#kable_styling("striped", "hover", full_width = F) %>%
#collapse_rows(columns = 1:2, latex_hline = "full") %>%
pass-> sampleBreakdown.html

sampleBreakdown.html

sampleBreakdown.html %>% cat(file = 'tables/sampleBreakdown.html')

Variability within vs between participants

psN2A.wuf <- phyloseq::distance(psN2A, method = "wunifrac")

psN2A %>% sample_data %>% .[1:5, 1:10]

##           pub_id    sex muVisit muVisitRank rx_code ct
## Sample-57    282   Male       9           1      T1  T
## Sample-58    282   Male      12           2      T1  T
## Sample-60    176 Female       9           1      T1  T
## Sample-61    176 Female      12           2      T1  T
## Sample-62    123 Female      12           1      T4  T
##           CD4._ANY.ENV.PTEG_Month_0 CD4._ANY.ENV.PTEG_Month_12
## Sample-57                     0.025                 0.02500000
## Sample-58                     0.025                 0.02500000
## Sample-60                     0.025                 0.06841053
## Sample-61                     0.025                 0.06841053
## Sample-62                     0.025                 0.42710324
##           CD4._ANY.ENV.PTEG_Month_6.5 IgA_gp41_Month_0
## Sample-57                  0.06146979            109.5
## Sample-58                  0.06146979            109.5
## Sample-60                  0.16678641            149.5
## Sample-61                  0.16678641            149.5
## Sample-62                  0.70855159             33.5

psN2A %>% sample_data %>% as("data.frame") %>% rownames_to_column(var = "Sample") %>% dplyr::select(Sample, pub_id) %>%
pass -> S2P

All.equal <- Vectorize(function(x,y){x == y})

## Convert distance matrix into long form matrix
psN2A.wuf %>% as.matrix %>% as.data.frame %>%
rownames_to_column(var = "SampleX")%>%
gather(key = "SampleY", value = "WufDist", -SampleX) %>%
left_join(S2P, by = c("SampleX" = "Sample")) %>% rename(pub_id_x = pub_id) %>%
left_join(S2P, by = c("SampleY" = "Sample")) %>% rename(pub_id_y = pub_id) %>%
mutate(SampleX = as.numeric(str_extract(SampleX, "[0-9][0-9]"))) %>%
mutate(SampleY = as.numeric(str_extract(SampleY, "[0-9][0-9]"))) %>%
## discard diagonal discard and upper half of triangular matrix
filter(SampleX < SampleY) %>%
mutate(isSamePerson = All.equal(pub_id_x, pub_id_y)) %>%
# ## discard cases where pub_id is unknown
# filter(is.finite(pub_id_x) & is.finite(pub_id_y))
pass -> AllWufDist

AllWufDist %>% head

##   SampleX SampleY   WufDist pub_id_x pub_id_y isSamePerson
## 1      57      58 0.3769328      282      282         TRUE
## 2      57      60 0.2754723      282      176        FALSE
## 3      58      60 0.3852072      282      176        FALSE
## 4      57      61 0.4105227      282      176        FALSE
## 5      58      61 0.3384203      282      176        FALSE
## 6      60      61 0.3868663      176      176         TRUE

#AllWufDist %>% ggplot(aes(x = isSamePerson, y = WufDist)) + geom_violin() + geom_dotplot(binaxis = "y", stackdir = "center", binwidth = .005)

Get Mean values for between participant and within participant weighted unifrac distances.

WufMeans <- AllWufDist %>% group_by(isSamePerson) %>% summarize(mean = mean(WufDist))
WufMeans

## # A tibble: 2 x 2
##   isSamePerson  mean
##   <lgl>        <dbl>
## 1 FALSE        0.497
## 2 TRUE         0.411

#SameAndDiff <- data.frame(comparison = c("different", "same"), WufMeans$mean)
SameAndDiff <- WufMeans %>% spread(key = isSamePerson, value = mean) %>% rename(between = 'FALSE', within = 'TRUE') %>% mutate(diff = between-within)
SameAndDiff

## # A tibble: 1 x 3
##   between within   diff
##     <dbl>  <dbl>  <dbl>
## 1   0.497  0.411 0.0860

#SameAndDiff[1,2] - SameAndDiff[2,2] # Difference in weighted unifrac dissimilarity between same and different partcipants

# Split the data
SamePersonWufDist <- AllWufDist %>% filter(isSamePerson) #%>% dplyr::select(WufDist)
DifferentPersonWufDist <- AllWufDist %>% filter(!isSamePerson)
set.seed(334)

bootsSame <- rsample::bootstraps(SamePersonWufDist, times = 10000)
bootsDifferent <- rsample::bootstraps(DifferentPersonWufDist, times = 10000)

mean_of_bootstrap <- function(split){
    locVals <- rsample::analysis(split)$WufDist
    mean(locVals)
}

boot_meansSame <- bootsSame %>% mutate(mean = map_dbl(splits, mean_of_bootstrap)) %>% dplyr::select(mean)

boot_meansDifferent <- bootsDifferent %>% mutate(mean = map_dbl(splits, mean_of_bootstrap)) %>% dplyr::select(mean)

boot_means <- bind_cols(within = boot_meansSame$mean, between = boot_meansDifferent$mean) %>%
mutate(isDifferentBigger = between>within, 
      DifMinusSame = between - within)
#boot_means

1- sum(boot_means$isDifferentBigger)/length(boot_means$isDifferentBigger)

## [1] 0.012

boot_means %>% ggplot(aes(x = DifMinusSame)) + geom_histogram()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

quantile(boot_means$DifMinusSame, c(0.025, 0.5, 0.975))

##       2.5%        50%      97.5% 
## 0.01273644 0.08721308 0.15521066

PermWufDist <- AllWufDist %>% modelr::permute(10000, WufDist)
mean_of_is_same <- function(df){df %>% as.data.frame %>% group_by(isSamePerson) %>% summarize(mean(WufDist)) %>% spread(isSamePerson, `mean(WufDist)`)}

test <- PermWufDist %>% pull(perm) %>% .[[1]] %>% as.data.frame
test %>% head

##   SampleX SampleY   WufDist pub_id_x pub_id_y isSamePerson
## 1      57      58 0.2863799      282      282         TRUE
## 2      57      60 0.5667858      282      176        FALSE
## 3      58      60 0.6974452      282      176        FALSE
## 4      57      61 0.8762536      282      176        FALSE
## 5      58      61 0.3909422      282      176        FALSE
## 6      60      61 0.8243214      176      176         TRUE

mean_of_is_same(test)

## # A tibble: 1 x 2
##   `FALSE` `TRUE`
##     <dbl>  <dbl>
## 1   0.496  0.462

SameAndDiff

## # A tibble: 1 x 3
##   between within   diff
##     <dbl>  <dbl>  <dbl>
## 1   0.497  0.411 0.0860

permutedMeans <- PermWufDist %>% mutate(meanvals = map(perm, mean_of_is_same)) %>% dplyr::select(meanvals) %>% unnest %>% 
rename(between = `FALSE`, within = `TRUE`) %>%
mutate(diff = between - within, absdif = abs(diff)) %>%
mutate(isExtreme = absdif >= SameAndDiff$diff, isExtreme1Tail = diff >= SameAndDiff$diff)
#permutedMeans %>% ggplot(aes(x = diff)) + geom_histogram() + geom_vline(xintercept = SameAndDiff$diff)

permutedMeans %>% ggplot(aes(x = diff)) + geom_histogram() + geom_vline(xintercept = SameAndDiff$diff)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#permutedMeans %>% mutate(isExtreme = absdif >= SameAndDiff$diff)
sum(permutedMeans$isExtreme) / length(permutedMeans$isExtreme)

## [1] 0.0153

sum(permutedMeans$isExtreme1Tail) / length(permutedMeans$isExtreme1Tail)

## [1] 0.0073

options(repr.plot.width=6, repr.plot.height= 6)
boot_means  %>% dplyr::select(within, between) %>% gather(key = "comparison", value = "WufDist") %>% ggplot(aes(x = comparison, y = WufDist)) + 
geom_dotplot(data = AllWufDist %>% mutate(isSamePerson2 = if_else(isSamePerson, "within", "between"))
             , aes(x = isSamePerson2, y = WufDist), binaxis = "y", stackdir = "center", binwidth = .01, colour = "gray40", fill = "white") + 
geom_violin(fill = NA) + 
geom_point(data = data.frame(comparison = c("within", "between"), WufDist = WufMeans$mean), aes(x = comparison, y = c(WufDist[2], WufDist[1])), shape = 22, fill = "black", size = 2)+
theme_bw() +
scale_x_discrete(limits = c("within", "between")) + labs(y = "Weighted Unifrac Distance") #

ggsave('figures/BetweenVsWithin.png')  

## Saving 7 x 5 in image

Variability at earliest sampling

psN2.wuf <- phyloseq::distance(psN2, method = "wunifrac")

psN2.pcoa <- capscale(psN2.wuf ~ 1)

# How much variance si explained by each weighted unifrac axis
# Note, ten axes cover 95% of the variance. 
# I'm not going to look beyond that for any test.
data.frame(eig = psN2.pcoa$CA$eig) %>%
rownames_to_column('axis') %>%
mutate(proportion = eig/sum(eig)) %>%
mutate(cumulative = cumsum(proportion))

##     axis         eig   proportion cumulative
## 1   MDS1 0.825438822 0.2888040980  0.2888041
## 2   MDS2 0.495612315 0.1734045743  0.4622087
## 3   MDS3 0.364619128 0.1275727475  0.5897814
## 4   MDS4 0.268169168 0.0938268864  0.6836083
## 5   MDS5 0.261743609 0.0915787154  0.7751870
## 6   MDS6 0.162405998 0.0568225245  0.8320095
## 7   MDS7 0.146142884 0.0511323948  0.8831419
## 8   MDS8 0.097842330 0.0342330228  0.9173750
## 9   MDS9 0.061705028 0.0215893226  0.9389643
## 10 MDS10 0.050336955 0.0176118673  0.9565762
## 11 MDS11 0.041072011 0.0143702533  0.9709464
## 12 MDS12 0.022613493 0.0079119971  0.9788584
## 13 MDS13 0.017745217 0.0062086872  0.9850671
## 14 MDS14 0.016323022 0.0057110902  0.9907782
## 15 MDS15 0.014802146 0.0051789668  0.9959571
## 16 MDS16 0.009677461 0.0033859449  0.9993431
## 17 MDS17 0.001877523 0.0006569069  1.0000000

my_breaks = c(1, 75, 250, 500, 1000,2000)
psN2 %>% mutate_phyloseq_sample(
                               mc41 = factor(medcode_hl(IgG_gp41_Month_0)),
                                                log120 = (IgG_Con.6.gp120.B_Month_12)) -> psN2_mod
psN2_mod%>%
sample_data() %>%
ggplot(aes(x = MDS1, y = MDS2)) + geom_point(aes(fill = mc41), size = 5, stroke = 1, shape = 21) +
coord_fixed(sqrt(psN2.pcoa$CA$eig[2]/psN2.pcoa$CA$eig[1])) +
viridis::scale_fill_viridis(name = 'gp41 Baseline', direction = -1, discrete = TRUE) +
#scale_colour_manual(name = 'gp41 Primary', values = c('black', 'grey70')) + 
theme_bw() -> wuford_gp41


psN2_mod %>%
sample_data() %>%
ggplot(aes(x = MDS1, y = MDS2)) + geom_point(aes(fill = log120), size = 5, stroke = 1, shape = 21) +
coord_fixed(sqrt(psN2.pcoa$CA$eig[2]/psN2.pcoa$CA$eig[1])) +
viridis::scale_fill_viridis(name = 'gp120 Final', direction = 1, trans = "sqrt",
                           breaks = my_breaks, labels = my_breaks) +
#scale_colour_manual(name = 'gp41 Primary', values = c('black', 'grey70')) + 
theme_bw() -> wuford_gp120

par <- options()
options(repr.plot.width=11, repr.plot.height= 4)
#g <- grid.arrange(wuford_gp41, wuford_gp120, ncol = 2)
g <- cowplot::plot_grid(wuford_gp41, wuford_gp120, ncol = 2, labels = c("A", "B"), label_size = 24)
g

#ggsave('figures/wunifrac_Agp41_Bgp120_pcoa.png', g, width = 8, height = 4)
#ggsave('figures/wunifrac_Agp41_Bgp120_pcoa.svg', g, width = 8, height = 4)
cowplot::save_plot('figures/wunifrac_Agp41_Bgp120_pcoa.png', g, base_width = 8, base_height = 4)
cowplot::save_plot('figures/wunifrac_Agp41_Bgp120_pcoa.svg', g, base_width = 8, base_height = 4)

Table 1

# export conditionally formatted table as html

colNames1 = c(' ' = 3, 'Kernel' = 2, 'MDS' = 4)
colNames2 = c('Type', 'Antigen', 'Month', 'P', 'Q', 'P', 'Q', 'R2', 'Coef' )

concisePermKernTable %>%
mutate(
    # this row needs to happen first, since the reformatting of the nother numbers makes them harder to call
    MDS1_Coef = cell_spec(format(MDS1_Coef, digits = 3), "html",
                           bold = ifelse(Kernel_P < 0.05, 
                                         T,
                                         F),
                          italic = ifelse(Kernel_P < 0.05 & MDS1_Coef < 0,
                                        T,
                                         F),
                          background = ifelse(Kernel_P < 0.05, ifelse(MDS1_Coef < 0, "lightsalmon", "lightblue"), "")
                         ),
    Kernel_P = cell_spec(format_round(Kernel_P, 3), "html",
                                  bold = ifelse(Kernel_P < 0.05, T, F),
                                  background = ifelse(Kernel_P < 0.05, 'yellow', '')
                                 ),
    Kernel_Q = cell_spec(format_round(Kernel_Q, 3), "html",
                                  bold = ifelse(Kernel_Q < 0.2, T, F),
                                  background = ifelse(Kernel_Q < 0.2, 'lightyellow', '')
                                 ),
     MDS1_P  = cell_spec(format_round(MDS1_P, 3), "html",
                                  bold = ifelse(MDS1_P < 0.05, T, F),
                                  background = ifelse(MDS1_P < 0.05, 'yellow', '')
                                 ),
    MDS1_Q = cell_spec(format_round(MDS1_Q, 3), "html",
                                  bold = ifelse(MDS1_Q < 0.2, T, F),
                                  background = ifelse(MDS1_Q < 0.2, 'lightyellow', '')
                                 ),
    #Month = cell_spec(format_round(Month,0), "html")
    Month = cell_spec(Month, "html")

    
      ) %>%
mutate(Antigen = gsub('ANY.ENV.PTEG', 'Any ENV PTEG', Antigen)) %>%
mutate(Antigen = gsub('gp70_B.CaseA_V1_V2', 'gp70 B.CaseA V1-V2', Antigen)) -> toTable

toTable %>%

kable("html", escape = F, digits = 3, align = 'c', col.names = colNames2) %>%
kable_styling("striped", "hover", full_width = F) %>%
add_header_above(colNames1) %>%
collapse_rows(columns = 1:2, latex_hline = "full") -> concisePermKernTable.html

concisePermKernTable.html

concisePermKernTable.html %>% cat(file = 'tables/concisePermkernTable.html')

Latex version of the same table

docHead <- "\\documentclass[12pt]{article} % use larger type; default would be 10pt

\\usepackage[utf8]{inputenc} % set input encoding (not needed with XeLaTeX)
\\usepackage{booktabs}
\\usepackage{longtable}
\\usepackage{array}
\\usepackage{multirow}
\\usepackage[table]{xcolor}
\\usepackage{wrapfig}
\\usepackage{float}
\\usepackage{colortbl}
\\usepackage{pdflscape}
\\usepackage{tabu}
\\usepackage{threeparttable}
\\usepackage{threeparttablex}
\\usepackage[normalem]{ulem}
\\usepackage{makecell}

\\definecolor{green}{rgb}{1, 1, .9}

\\begin{document}
"

docTail <- "\\end{document}
"

# Make latex table

concisePermKernTable %>%
mutate(
    # this row needs to happen first, since the reformatting of the nother numbers makes them harder to call
    MDS1_Coef = cell_spec(format(MDS1_Coef, digits = 3), "latex",
                           bold = ifelse(MDS1_P < 0.05, 
                                         T,
                                         F),
                          italic = ifelse(MDS1_P < 0.05 & MDS1_Coef < 0,
                                        T,
                                         F)),
    Kernel_P = cell_spec(format_round(Kernel_P, 3), "latex",
                                  bold = ifelse(Kernel_P < 0.05, T, F),
                                  background = ifelse(Kernel_P < 0.05, 'yellow', 'white')
                                 ),
    Kernel_Q = cell_spec(format_round(Kernel_Q, 3), "latex",
                                  bold = ifelse(Kernel_Q < 0.2, T, F),
                                  background = ifelse(Kernel_Q < 0.2, 'green', 'white')
                                 ),
     MDS1_P  = cell_spec(format_round(MDS1_P, 3), "latex",
                                  bold = ifelse(MDS1_P < 0.05, T, F),
                                  background = ifelse(MDS1_P < 0.05, 'yellow', 'white')
                                 ),
    MDS1_Q = cell_spec(format_round(MDS1_Q, 3), "latex",
                                  bold = ifelse(MDS1_Q < 0.2, T, F),
                                  background = ifelse(MDS1_Q < 0.2, 'green', 'white')
                                 ),
    #Month = cell_spec(format_round(Month,0), "html")
    Month = cell_spec(Month, "latex")

    
      ) %>%
mutate(Antigen = gsub('ANY.ENV.PTEG', 'Any ENV PTEG', Antigen)) %>%
mutate(Antigen = gsub('gp70_B.CaseA_V1_V2', 'gp70 B.CaseA V1-V2', Antigen)) -> toTable

toTable %>% 
kable("latex", escape = F, digits = 3, align = 'c', col.names = colNames2, booktabs = T) %>%
kable_styling(position = "left") %>%

add_header_above(colNames1) %>%
collapse_rows(columns = 1:2, latex_hline = "full") %>%
pass -> concisePermKernTable.tex

# Print latex table to tex file

cat(docHead, concisePermKernTable.tex, docTail, file = 'tables/concisePermkernTable1.tex')

concisePermKernTable %>% filter(Kernel_P < 0.05) -> shortPermkernTable
shortPermkernTable

##   Type            Antigen Month Kernel_P Kernel_Q      MDS1_P     MDS1_Q
## 1  IgG      Con.6.gp120.B   6.5  0.00431   0.0862 0.001759476 0.01064944
## 2  IgG      Con.6.gp120.B  12.0  0.03584   0.1538 0.010022385 0.01740207
## 3  IgG               gp41   0.0  0.04614   0.1538 0.014375638 0.01740207
## 4  IgG               gp41   6.5  0.04457   0.1538 0.030461387 0.03072856
## 5  IgG gp70_B.CaseA_V1_V2  12.0  0.03465   0.1538 0.014277612 0.01740207
## 6  IgG         ZM96.gp140   6.5  0.01671   0.1538 0.008971440 0.01740207
##   GlmMDS1_R2 MDS1_Coef
## 1  0.4969907 -3.103554
## 2  0.3612921 -2.289438
## 3  0.3312048  2.184936
## 4  0.2667186 -1.808614
## 5  0.3317812 -2.134684
## 6  0.3704046 -2.338385

write.csv(format(shortPermkernTable, digits = 3), 'tables/shortPermkernTable.csv')

Q-Q Plot

Of kernel regression P values

my_runif <- function(Len){
    loc_runif <- runif(n = Len)
    sort_loc_runif <- sort(loc_runif)
    data.frame(case = 1:Len, relement = sort_loc_runif)
}

calc_bound <- function(df, bound){
    quantile(df$relement, bound)
}

make_qqdata <- function(pvec, nboot = 10000){
    locP <- pvec

LP <- length(locP)
#1:10 %>% map(runif, n = LP)
sortedP <- sort(locP)
exp2 <- 1:length(locP)/length(locP)
sortedExpP <- sort(exp2)

random_pvalues <- data.frame(iter = 1:nboot) %>%
mutate(rand = map(iter, ~my_runif(Len = LP))) %>% unnest %>%
nest(-case) %>%
mutate(lb = map_dbl(data, ~calc_bound(., 0.025)),
      ub = map_dbl(data, ~calc_bound(.,0.975))) %>% dplyr::select (-data)

qqdata <- bind_cols(sortedP = sortedP, sortedExpP = sortedExpP, random_pvalues)

return(qqdata)
}

qqdata_permKernMir <- make_qqdata(permKernTable$mir.P)

qqdata_permKernMir %>% str

## Classes 'tbl_df', 'tbl' and 'data.frame':    20 obs. of  5 variables:
##  $ sortedP   : num  0.00431 0.01671 0.03465 0.03584 0.04457 ...
##  $ sortedExpP: num  0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 ...
##  $ case      : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ lb        : num  0.00123 0.01314 0.03293 0.05827 0.08676 ...
##  $ ub        : num  0.169 0.247 0.317 0.38 0.441 ...

options(repr.plot.width=6, repr.plot.height= 6)
qqdata_permKernMir %>% ggplot(aes(x = sortedExpP, y = sortedP)) + geom_point(shape = 1) + geom_abline(intercept = 0, slope = 1) + 
labs(x = expression(paste("Expected ", italic("p"),"-value")), y = expression(paste("Observed ", italic("p"),"-value")))

ggsave('figures/qq_permKernMir.png')

## Saving 7 x 5 in image

qqdata_permKernMir %>% ggplot(aes(x = sortedExpP, y = sortedP)) + geom_point(shape = 1) + geom_abline(intercept = 0, slope = 1) + geom_line(aes(y = lb), colour = "grey40") + geom_line(aes(y = ub), colour = "grey40") + scale_y_log10() + scale_x_log10() + 
  labs(x = expression(paste("Expected ", italic("p"),"-value")), y = expression(paste("Observed ", italic("p"),"-value")))

We see from the qqplot that we generally fall below the 1:1 line, suggesting that our P values are lower than we would expect from chance. To make the grey lines, I sampled 95% confidence intervals of random draws from the uniform distribution (essentially random p values). The data points sometimes, but not always fall below the lower one of these. However, we care about all of the p values, not wheter eaach individual is more unusual than we would expect by chance. I’ll replort just the regular qqplot as a supplemental figure.

As above, but this time with gaussian - continuous dependent variables

ptm = proc.time()
tps <- CapVar(test.immune1, nperm = 9999, transformation = function(x){jac_box_cox_2(x)}, family = 'gaussian')
proc.time() - ptm

##    user  system elapsed 
##   0.648   0.000   0.648

tps

##      caps.P adonisP  mir.P   caps.F    caps.R2     wuf1.P   wuf1.DR
## MDS1 0.2222  0.2538 0.2194 1.311871 0.06520801 0.03795074 0.1848023
##      wuf1.McFadden accuracy wuf1.coef
## MDS1     0.1969076        0  0.700859

# Run above function against every relevant variable.
ptm <- proc.time()

use.immune %>%
filter(ct == 'T') %>%
group_by(type, antigen, month) %>%
do(data.frame(CapVar(., nperm = jnperm,
                     transformation = function(x){jac_box_cox_2(x)},
                     family = 'gaussian'))) -> permKernTableGaus
permKernTableGaus

## # A tibble: 20 x 13
## # Groups:   type, antigen, month [20]
##    type  antigen month caps.P adonisP   mir.P caps.F caps.R2  wuf1.P
##    <chr> <chr>   <dbl>  <dbl>   <dbl>   <dbl>  <dbl>   <dbl>   <dbl>
##  1 CD4+  ANY.EN…   6.5 0.603  6.02e-1 0.602    0.827  0.0421 5.97e-1
##  2 CD4+  ANY.EN…  12   0.185  1.94e-1 0.185    1.36   0.0707 3.29e-1
##  3 IgA   gp41      0   0.990  9.88e-1 0.989    0.311  0.0163 6.33e-1
##  4 IgA   gp41      6.5 0.411  4.07e-1 0.411    1.00   0.0505 6.01e-1
##  5 IgA   gp41     12   0.655  6.52e-1 0.651    0.755  0.0406 5.04e-1
##  6 IgA   p24       0   0.953  9.53e-1 0.954    0.417  0.0217 4.33e-1
##  7 IgA   p24       6.5 0.976  9.73e-1 0.976    0.373  0.0194 8.57e-1
##  8 IgA   p24      12   0.574  5.75e-1 0.575    0.852  0.0456 3.09e-1
##  9 IgG   Con.6.…   6.5 0.0736 7.60e-2 0.0731   1.72   0.0839 3.16e-2
## 10 IgG   Con.6.…  12   0.0007 6.00e-4 0.00068  3.55   0.159  1.43e-5
## 11 IgG   gp41      0   0.221  2.48e-1 0.221    1.31   0.0652 3.80e-2
## 12 IgG   gp41      6.5 0.220  2.25e-1 0.219    1.29   0.0643 7.98e-2
## 13 IgG   gp41     12   0.266  2.71e-1 0.263    1.23   0.0613 1.73e-1
## 14 IgG   gp70_B…   6.5 0.287  2.93e-1 0.282    1.19   0.0595 2.06e-1
## 15 IgG   gp70_B…  12   0.110  1.12e-1 0.111    1.55   0.0763 8.34e-2
## 16 IgG   p24       0   0.754  7.45e-1 0.753    0.636  0.0327 9.07e-1
## 17 IgG   p24       6.5 0.0268 2.61e-2 0.0272   2.15   0.108  5.30e-2
## 18 IgG   p24      12   0.411  4.14e-1 0.413    1.02   0.0513 1.33e-1
## 19 IgG   ZM96.g…   6.5 0.0293 2.96e-2 0.0290   2.24   0.106  1.69e-2
## 20 IgG   ZM96.g…  12   0.204  2.00e-1 0.198    1.32   0.0657 7.77e-2
## # ... with 4 more variables: wuf1.DR <dbl>, wuf1.McFadden <dbl>,
## #   accuracy <dbl>, wuf1.coef <dbl>

proc.time() - ptm

##    user  system elapsed 
## 117.617   0.532 118.160

# Clean up so we just see the results of the kernel regression 
concisePermKernTableGaus <- permKernTableGaus %>% ungroup %>%
mutate(Kernel_Q = p2q(mir.P), MDS1_Q = p2q(wuf1.P)) %>%
dplyr::select(Type = type, Antigen = antigen, Month = month, Kernel_P = mir.P, Kernel_Q,
              MDS1_P = wuf1.P, MDS1_Q, MDS1_R2 = wuf1.McFadden, MDS1_Coef = wuf1.coef) %>%
as.data.frame %>% 
pass 

concisePermKernTableGaus

##    Type            Antigen Month Kernel_P  Kernel_Q       MDS1_P
## 1  CD4+       ANY.ENV.PTEG   6.5  0.60213 0.8028400 5.971764e-01
## 2  CD4+       ANY.ENV.PTEG  12.0  0.18544 0.4920667 3.291117e-01
## 3   IgA               gp41   0.0  0.98896 0.9889600 6.330939e-01
## 4   IgA               gp41   6.5  0.41053 0.6358769 6.006684e-01
## 5   IgA               gp41  12.0  0.65053 0.8131625 5.043440e-01
## 6   IgA                p24   0.0  0.95380 0.9889600 4.331569e-01
## 7   IgA                p24   6.5  0.97610 0.9889600 8.573830e-01
## 8   IgA                p24  12.0  0.57478 0.8028400 3.091250e-01
## 9   IgG      Con.6.gp120.B   6.5  0.07311 0.3655500 3.164215e-02
## 10  IgG      Con.6.gp120.B  12.0  0.00068 0.0136000 1.431333e-05
## 11  IgG               gp41   0.0  0.22143 0.4920667 3.795074e-02
## 12  IgG               gp41   6.5  0.21873 0.4920667 7.975175e-02
## 13  IgG               gp41  12.0  0.26339 0.5136000 1.731298e-01
## 14  IgG gp70_B.CaseA_V1_V2   6.5  0.28248 0.5136000 2.062894e-01
## 15  IgG gp70_B.CaseA_V1_V2  12.0  0.11097 0.4438800 8.341959e-02
## 16  IgG                p24   0.0  0.75279 0.8856353 9.069707e-01
## 17  IgG                p24   6.5  0.02725 0.1932667 5.296858e-02
## 18  IgG                p24  12.0  0.41332 0.6358769 1.334352e-01
## 19  IgG         ZM96.gp140   6.5  0.02899 0.1932667 1.691120e-02
## 20  IgG         ZM96.gp140  12.0  0.19822 0.4920667 7.767448e-02
##          MDS1_Q      MDS1_R2   MDS1_Coef
## 1  0.3441069304 0.0154346775 -0.19622223
## 2  0.2476839187 0.0535522002 -0.35622365
## 3  0.3441069304 0.0126283362  0.17748943
## 4  0.3441069304 0.0151466567 -0.19438279
## 5  0.3289526950 0.0257545280  0.25159262
## 6  0.3027018097 0.0333728299  0.28853330
## 7  0.4414881726 0.0018079718 -0.06715761
## 8  0.2476839187 0.0579135356 -0.37727770
## 9  0.0928235317 0.2083288787 -0.72089857
## 10 0.0001400356 0.5303144044 -1.15018086
## 11 0.0928235317 0.1969076139  0.70085904
## 12 0.1020177828 0.1482127646 -0.60805411
## 13 0.1693829682 0.0948033344 -0.48630779
## 14 0.1834771570 0.0826277473 -0.45400682
## 15 0.1020177828 0.1451699474 -0.60178005
## 16 0.4436710261 0.0007652868  0.04369297
## 17 0.1020177828 0.1835313650 -0.78663031
## 18 0.1450526857 0.1129112051 -0.53072304
## 19 0.0827260298 0.2460702885 -0.78348198
## 20 0.1020177828 0.1499943359 -0.61169771

write.csv(format(concisePermKernTableGaus, digits = 3), 'tables/concisePermkernTableGaus.csv')

# export conditionally formatted table as html
concisePermKernTableGaus %>%
mutate(
    # this row needs to happen first, since the reformatting of the nother numbers makes them harder to call
       MDS1_Coef = cell_spec(format_round(MDS1_Coef, digits = 3), "html",
                           bold = ifelse(Kernel_P < 0.05, 
                                         T,
                                         F),
                          italic = ifelse(Kernel_P < 0.05 & MDS1_Coef < 0,
                                        T,
                                         F),
                          background = ifelse(Kernel_P < 0.05, ifelse(MDS1_Coef < 0, "lightsalmon", "lightblue"), "")
                         ),
    Kernel_P = cell_spec(format_round(Kernel_P, 3), "html",
                                  bold = ifelse(Kernel_P < 0.05, T, F),
                                  background = ifelse(Kernel_P < 0.05, 'yellow', '')
                                 ),
    Kernel_Q = cell_spec(format_round(Kernel_Q, 3), "html",
                                  bold = ifelse(Kernel_Q < 0.2, T, F),
                                  background = ifelse(Kernel_Q < 0.2, 'lightyellow', '')
                                 ),
     MDS1_P  = cell_spec(format_round(MDS1_P, 3), "html",
                                  bold = ifelse(MDS1_P < 0.05, T, F),
                                  background = ifelse(MDS1_P < 0.05, 'yellow', '')
                                 ),
    MDS1_Q = cell_spec(format_round(MDS1_Q, 3), "html",
                                  bold = ifelse(MDS1_Q < 0.2, T, F),
                                  background = ifelse(MDS1_Q < 0.2, 'lightyellow', '')
                                 ),
    Month = cell_spec(Month, "html")

    
      ) %>%

mutate(Antigen = gsub('ANY.ENV.PTEG', 'Any ENV PTEG', Antigen)) %>%
mutate(Antigen = gsub('gp70_B.CaseA_V1_V2', 'gp70 B.CaseA V1-V2', Antigen)) -> toTable

toTable %>%

kable("html", escape = F, digits = 3, align = 'c', col.names = colNames2) %>%
kable_styling("striped", "hover", full_width = F) %>%
add_header_above(colNames1) %>%
collapse_rows(columns = 1:2, latex_hline = "full") -> concisePermKernTableGaus.html

concisePermKernTableGaus.html

concisePermKernTableGaus.html %>% cat(file = 'tables/concisePermkernTableGaus.html')

concisePermKernTableGaus %>% filter(Kernel_P < 0.05) -> shortPermkernTableGaus
shortPermkernTableGaus

##   Type       Antigen Month Kernel_P  Kernel_Q       MDS1_P       MDS1_Q
## 1  IgG Con.6.gp120.B  12.0  0.00068 0.0136000 1.431333e-05 0.0001400356
## 2  IgG           p24   6.5  0.02725 0.1932667 5.296858e-02 0.1020177828
## 3  IgG    ZM96.gp140   6.5  0.02899 0.1932667 1.691120e-02 0.0827260298
##     MDS1_R2  MDS1_Coef
## 1 0.5303144 -1.1501809
## 2 0.1835314 -0.7866303
## 3 0.2460703 -0.7834820

write.csv(format(shortPermkernTableGaus, digits = 3), 'tables/shortPermkernTable.csv')

# Make latex table

concisePermKernTableGaus %>%
mutate(
    # this row needs to happen first, since the reformatting of the nother numbers makes them harder to call
    MDS1_Coef = cell_spec(format(MDS1_Coef, digits = 3), "html",
                           bold = ifelse(Kernel_P < 0.05, 
                                         T,
                                         F),
                          italic = ifelse(Kernel_P < 0.05 & MDS1_Coef < 0,
                                        T,
                                         F),
                          background = ifelse(Kernel_P < 0.05, ifelse(MDS1_Coef < 0, "lightsalmon", "lightblue"), "")
                         ),
    Kernel_P = cell_spec(format_round(Kernel_P, 3), "latex",
                                  bold = ifelse(Kernel_P < 0.05, T, F),
                                  background = ifelse(Kernel_P < 0.05, 'yellow', 'white')
                                 ),
    Kernel_Q = cell_spec(format_round(Kernel_Q, 3), "latex",
                                  bold = ifelse(Kernel_Q < 0.2, T, F),
                                  background = ifelse(Kernel_Q < 0.2, 'green', 'white')
                                 ),
     MDS1_P  = cell_spec(format_round(MDS1_P, 3), "latex",
                                  bold = ifelse(MDS1_P < 0.05, T, F),
                                  background = ifelse(MDS1_P < 0.05, 'yellow', 'white')
                                 ),
    MDS1_Q = cell_spec(format_round(MDS1_Q, 3), "latex",
                                  bold = ifelse(MDS1_Q < 0.2, T, F),
                                  background = ifelse(MDS1_Q < 0.2, 'green', 'white')
                                 ),
    #Month = cell_spec(format_round(Month,0), "html")
    Month = cell_spec(Month, "latex")

    
      ) %>%
mutate(Antigen = gsub('ANY.ENV.PTEG', 'Any ENV PTEG', Antigen)) %>%
mutate(Antigen = gsub('gp70_B.CaseA_V1_V2', 'gp70 B.CaseA V1-V2', Antigen)) -> toTable

toTable %>% 
kable("latex", escape = F, digits = 3, align = 'c', col.names = colNames2, booktabs = T) %>%
kable_styling(position = "left") %>%

add_header_above(colNames1) %>%
collapse_rows(columns = 1:2, latex_hline = "full") %>%
pass -> concisePermKernTableGaus.tex

# Print latex table to tex file

cat(docHead, concisePermKernTableGaus.tex, docTail, file = 'tables/concisePermkernTableGaus.tex')

### QQplot for kernel regression data

qqdata_permKernMirGaus <- make_qqdata(permKernTableGaus$mir.P)

options(repr.plot.width=6, repr.plot.height= 6)
qqdata_permKernMirGaus %>% ggplot(aes(x = sortedExpP, y = sortedP)) + geom_point(shape = 1) + geom_abline(intercept = 0, slope = 1) + 
labs(x = expression(paste("Expected ", italic("p"),"-value")), y = expression(paste("Observed ", italic("p"),"-value")))

ggsave('figures/qq_permKernMirGaus.png')

## Saving 7 x 5 in image

Chi Squared test for statistical associations between each pair of immune variables

use.immune %>% dplyr::select(pub_id, visitno, type, antigen, mag) -> tmp
full_join(tmp, tmp, by = 'pub_id') %>% 
group_by(visitno.x, type.x, antigen.x, visitno.y, type.y, antigen.y) %>%
nest %>%
mutate(x2 = map(data, function(df){unwarn(chisq.test(df$mag.x, df$mag.y))})) %>%
mutate(glance = map(x2, glance)) %>%
dplyr::select(-data, -x2) %>%
unnest(glance) %>%
#mutate(q.value = p2q(p.value)) %>% # reurns NaNs
pass -> compareImmuneX2

compareImmuneX2 %>% filter(
    type.x == 'IgG' &
    antigen.x == 'gp41' &
    type.y == 'IgG' &
    antigen.y == 'gp41'
)

## # A tibble: 9 x 10
##   visitno.x type.x antigen.x visitno.y type.y antigen.y statistic p.value
##       <int> <chr>  <chr>         <int> <chr>  <chr>         <dbl>   <dbl>
## 1         9 IgG    gp41              9 IgG    gp41            420  0.236 
## 2         9 IgG    gp41             12 IgG    gp41            378  0.247 
## 3         9 IgG    gp41              2 IgG    gp41            420  0.236 
## 4        12 IgG    gp41              9 IgG    gp41            378  0.247 
## 5        12 IgG    gp41             12 IgG    gp41            378  0.0207
## 6        12 IgG    gp41              2 IgG    gp41            378  0.247 
## 7         2 IgG    gp41              9 IgG    gp41            420  0.236 
## 8         2 IgG    gp41             12 IgG    gp41            378  0.247 
## 9         2 IgG    gp41              2 IgG    gp41            420  0.236 
## # ... with 2 more variables: parameter <int>, method <chr>

compareImmuneX2 %>%
filter(type.x == 'IgG' & type.y == 'IgG' & antigen.x != antigen.y) %>%
write_csv('tables/chisq_IgG_comparasons.csv')

Table S2

allMDS %>%
kable("html", escape = F, digits = 3, align = 'c') %>%
collapse_rows(columns = 1:2, latex_hline = "full") %>%
as.character() -> allMDS.html

allMDS %>%
kable("html", escape = F, digits = 3, align = 'c') %>%
collapse_rows(columns = 1:2, latex_hline = "full")

allMDS.html

## [1] "<table>\n <thead>\n  <tr>\n   <th style=\"text-align:center;\"> Type </th>\n   <th style=\"text-align:center;\"> Antigen </th>\n   <th style=\"text-align:center;\"> Month </th>\n   <th style=\"text-align:center;\"> MDS1 </th>\n   <th style=\"text-align:center;\"> MDS2 </th>\n   <th style=\"text-align:center;\"> MDS3 </th>\n   <th style=\"text-align:center;\"> MDS4 </th>\n   <th style=\"text-align:center;\"> MDS5 </th>\n   <th style=\"text-align:center;\"> MDS6 </th>\n   <th style=\"text-align:center;\"> MDS7 </th>\n   <th style=\"text-align:center;\"> MDS8 </th>\n   <th style=\"text-align:center;\"> MDS9 </th>\n   <th style=\"text-align:center;\"> MDS10 </th>\n  </tr>\n </thead>\n<tbody>\n  <tr>\n   <td style=\"text-align:center;vertical-align: middle !important;\" rowspan=\"6\"> IgA </td>\n   <td style=\"text-align:center;vertical-align: middle !important;\" rowspan=\"3\"> gp41 </td>\n   <td style=\"text-align:center;\"> 0.0 </td>\n   <td style=\"text-align:center;\"> 0.653 </td>\n   <td style=\"text-align:center;\"> 0.607 </td>\n   <td style=\"text-align:center;\"> 0.401 </td>\n   <td style=\"text-align:center;\"> 0.702 </td>\n   <td style=\"text-align:center;\"> 0.786 </td>\n   <td style=\"text-align:center;\"> 0.650 </td>\n   <td style=\"text-align:center;\"> 0.242 </td>\n   <td style=\"text-align:center;\"> 0.349 </td>\n   <td style=\"text-align:center;\"> 0.038 </td>\n   <td style=\"text-align:center;\"> 0.700 </td>\n  </tr>\n  <tr>\n   \n   \n   <td style=\"text-align:center;\"> 6.5 </td>\n   <td style=\"text-align:center;\"> 0.168 </td>\n   <td style=\"text-align:center;\"> 0.368 </td>\n   <td style=\"text-align:center;\"> 0.701 </td>\n   <td style=\"text-align:center;\"> 0.954 </td>\n   <td style=\"text-align:center;\"> 0.112 </td>\n   <td style=\"text-align:center;\"> 0.142 </td>\n   <td style=\"text-align:center;\"> 0.456 </td>\n   <td style=\"text-align:center;\"> 0.475 </td>\n   <td style=\"text-align:center;\"> 0.789 </td>\n   <td style=\"text-align:center;\"> 0.864 </td>\n  </tr>\n  <tr>\n   \n   \n   <td style=\"text-align:center;\"> 12.0 </td>\n   <td style=\"text-align:center;\"> 0.855 </td>\n   <td style=\"text-align:center;\"> 0.201 </td>\n   <td style=\"text-align:center;\"> 0.462 </td>\n   <td style=\"text-align:center;\"> 0.408 </td>\n   <td style=\"text-align:center;\"> 0.994 </td>\n   <td style=\"text-align:center;\"> 0.798 </td>\n   <td style=\"text-align:center;\"> 0.469 </td>\n   <td style=\"text-align:center;\"> 0.106 </td>\n   <td style=\"text-align:center;\"> 0.079 </td>\n   <td style=\"text-align:center;\"> 0.232 </td>\n  </tr>\n  <tr>\n   \n   <td style=\"text-align:center;vertical-align: middle !important;\" rowspan=\"3\"> p24 </td>\n   <td style=\"text-align:center;\"> 0.0 </td>\n   <td style=\"text-align:center;\"> 0.335 </td>\n   <td style=\"text-align:center;\"> 0.629 </td>\n   <td style=\"text-align:center;\"> 0.520 </td>\n   <td style=\"text-align:center;\"> 0.939 </td>\n   <td style=\"text-align:center;\"> 0.666 </td>\n   <td style=\"text-align:center;\"> 0.657 </td>\n   <td style=\"text-align:center;\"> 0.148 </td>\n   <td style=\"text-align:center;\"> 0.853 </td>\n   <td style=\"text-align:center;\"> 0.921 </td>\n   <td style=\"text-align:center;\"> 0.517 </td>\n  </tr>\n  <tr>\n   \n   \n   <td style=\"text-align:center;\"> 6.5 </td>\n   <td style=\"text-align:center;\"> 0.652 </td>\n   <td style=\"text-align:center;\"> 0.281 </td>\n   <td style=\"text-align:center;\"> 0.681 </td>\n   <td style=\"text-align:center;\"> 0.414 </td>\n   <td style=\"text-align:center;\"> 0.438 </td>\n   <td style=\"text-align:center;\"> 0.599 </td>\n   <td style=\"text-align:center;\"> 0.999 </td>\n   <td style=\"text-align:center;\"> 0.199 </td>\n   <td style=\"text-align:center;\"> 0.881 </td>\n   <td style=\"text-align:center;\"> 0.351 </td>\n  </tr>\n  <tr>\n   \n   \n   <td style=\"text-align:center;\"> 12.0 </td>\n   <td style=\"text-align:center;\"> 0.398 </td>\n   <td style=\"text-align:center;\"> 0.724 </td>\n   <td style=\"text-align:center;\"> 0.227 </td>\n   <td style=\"text-align:center;\"> 0.202 </td>\n   <td style=\"text-align:center;\"> 0.128 </td>\n   <td style=\"text-align:center;\"> 0.465 </td>\n   <td style=\"text-align:center;\"> 0.445 </td>\n   <td style=\"text-align:center;\"> 0.156 </td>\n   <td style=\"text-align:center;\"> 0.516 </td>\n   <td style=\"text-align:center;\"> 0.159 </td>\n  </tr>\n  <tr>\n   <td style=\"text-align:center;vertical-align: middle !important;\" rowspan=\"12\"> IgG </td>\n   <td style=\"text-align:center;vertical-align: middle !important;\" rowspan=\"2\"> Con.6.gp120.B </td>\n   <td style=\"text-align:center;\"> 6.5 </td>\n   <td style=\"text-align:center;\"> 0.017 </td>\n   <td style=\"text-align:center;\"> 0.536 </td>\n   <td style=\"text-align:center;\"> 0.420 </td>\n   <td style=\"text-align:center;\"> 0.890 </td>\n   <td style=\"text-align:center;\"> 0.325 </td>\n   <td style=\"text-align:center;\"> 0.427 </td>\n   <td style=\"text-align:center;\"> 0.853 </td>\n   <td style=\"text-align:center;\"> 0.806 </td>\n   <td style=\"text-align:center;\"> 0.634 </td>\n   <td style=\"text-align:center;\"> 0.250 </td>\n  </tr>\n  <tr>\n   \n   \n   <td style=\"text-align:center;\"> 12.0 </td>\n   <td style=\"text-align:center;\"> 0.031 </td>\n   <td style=\"text-align:center;\"> 0.878 </td>\n   <td style=\"text-align:center;\"> 0.377 </td>\n   <td style=\"text-align:center;\"> 0.635 </td>\n   <td style=\"text-align:center;\"> 0.989 </td>\n   <td style=\"text-align:center;\"> 0.540 </td>\n   <td style=\"text-align:center;\"> 0.569 </td>\n   <td style=\"text-align:center;\"> 0.374 </td>\n   <td style=\"text-align:center;\"> 0.686 </td>\n   <td style=\"text-align:center;\"> 0.299 </td>\n  </tr>\n  <tr>\n   \n   <td style=\"text-align:center;vertical-align: middle !important;\" rowspan=\"3\"> gp41 </td>\n   <td style=\"text-align:center;\"> 0.0 </td>\n   <td style=\"text-align:center;\"> 0.040 </td>\n   <td style=\"text-align:center;\"> 0.894 </td>\n   <td style=\"text-align:center;\"> 0.304 </td>\n   <td style=\"text-align:center;\"> 0.967 </td>\n   <td style=\"text-align:center;\"> 0.442 </td>\n   <td style=\"text-align:center;\"> 0.231 </td>\n   <td style=\"text-align:center;\"> 0.504 </td>\n   <td style=\"text-align:center;\"> 0.595 </td>\n   <td style=\"text-align:center;\"> 0.856 </td>\n   <td style=\"text-align:center;\"> 0.509 </td>\n  </tr>\n  <tr>\n   \n   \n   <td style=\"text-align:center;\"> 6.5 </td>\n   <td style=\"text-align:center;\"> 0.057 </td>\n   <td style=\"text-align:center;\"> 0.274 </td>\n   <td style=\"text-align:center;\"> 0.226 </td>\n   <td style=\"text-align:center;\"> 0.891 </td>\n   <td style=\"text-align:center;\"> 0.858 </td>\n   <td style=\"text-align:center;\"> 0.165 </td>\n   <td style=\"text-align:center;\"> 0.601 </td>\n   <td style=\"text-align:center;\"> 0.426 </td>\n   <td style=\"text-align:center;\"> 0.212 </td>\n   <td style=\"text-align:center;\"> 0.606 </td>\n  </tr>\n  <tr>\n   \n   \n   <td style=\"text-align:center;\"> 12.0 </td>\n   <td style=\"text-align:center;\"> 0.779 </td>\n   <td style=\"text-align:center;\"> 0.310 </td>\n   <td style=\"text-align:center;\"> 0.443 </td>\n   <td style=\"text-align:center;\"> 0.228 </td>\n   <td style=\"text-align:center;\"> 0.680 </td>\n   <td style=\"text-align:center;\"> 0.268 </td>\n   <td style=\"text-align:center;\"> 0.995 </td>\n   <td style=\"text-align:center;\"> 0.039 </td>\n   <td style=\"text-align:center;\"> 0.357 </td>\n   <td style=\"text-align:center;\"> 0.054 </td>\n  </tr>\n  <tr>\n   \n   <td style=\"text-align:center;vertical-align: middle !important;\" rowspan=\"2\"> gp70 B.CaseA V1-V2 </td>\n   <td style=\"text-align:center;\"> 6.5 </td>\n   <td style=\"text-align:center;\"> 0.621 </td>\n   <td style=\"text-align:center;\"> 0.853 </td>\n   <td style=\"text-align:center;\"> 0.300 </td>\n   <td style=\"text-align:center;\"> 0.458 </td>\n   <td style=\"text-align:center;\"> 0.507 </td>\n   <td style=\"text-align:center;\"> 0.690 </td>\n   <td style=\"text-align:center;\"> 0.506 </td>\n   <td style=\"text-align:center;\"> 0.161 </td>\n   <td style=\"text-align:center;\"> 0.039 </td>\n   <td style=\"text-align:center;\"> 0.880 </td>\n  </tr>\n  <tr>\n   \n   \n   <td style=\"text-align:center;\"> 12.0 </td>\n   <td style=\"text-align:center;\"> 0.037 </td>\n   <td style=\"text-align:center;\"> 0.665 </td>\n   <td style=\"text-align:center;\"> 0.808 </td>\n   <td style=\"text-align:center;\"> 0.447 </td>\n   <td style=\"text-align:center;\"> 0.318 </td>\n   <td style=\"text-align:center;\"> 0.296 </td>\n   <td style=\"text-align:center;\"> 0.743 </td>\n   <td style=\"text-align:center;\"> 0.143 </td>\n   <td style=\"text-align:center;\"> 0.403 </td>\n   <td style=\"text-align:center;\"> 0.276 </td>\n  </tr>\n  <tr>\n   \n   <td style=\"text-align:center;vertical-align: middle !important;\" rowspan=\"3\"> p24 </td>\n   <td style=\"text-align:center;\"> 0.0 </td>\n   <td style=\"text-align:center;\"> 0.451 </td>\n   <td style=\"text-align:center;\"> 0.231 </td>\n   <td style=\"text-align:center;\"> 0.986 </td>\n   <td style=\"text-align:center;\"> 0.088 </td>\n   <td style=\"text-align:center;\"> 0.149 </td>\n   <td style=\"text-align:center;\"> 0.051 </td>\n   <td style=\"text-align:center;\"> 0.882 </td>\n   <td style=\"text-align:center;\"> 0.507 </td>\n   <td style=\"text-align:center;\"> 0.384 </td>\n   <td style=\"text-align:center;\"> 0.439 </td>\n  </tr>\n  <tr>\n   \n   \n   <td style=\"text-align:center;\"> 6.5 </td>\n   <td style=\"text-align:center;\"> 0.404 </td>\n   <td style=\"text-align:center;\"> 0.828 </td>\n   <td style=\"text-align:center;\"> 0.243 </td>\n   <td style=\"text-align:center;\"> 0.250 </td>\n   <td style=\"text-align:center;\"> 0.689 </td>\n   <td style=\"text-align:center;\"> 0.082 </td>\n   <td style=\"text-align:center;\"> 0.180 </td>\n   <td style=\"text-align:center;\"> 0.846 </td>\n   <td style=\"text-align:center;\"> 0.360 </td>\n   <td style=\"text-align:center;\"> 0.101 </td>\n  </tr>\n  <tr>\n   \n   \n   <td style=\"text-align:center;\"> 12.0 </td>\n   <td style=\"text-align:center;\"> 0.182 </td>\n   <td style=\"text-align:center;\"> 0.721 </td>\n   <td style=\"text-align:center;\"> 0.237 </td>\n   <td style=\"text-align:center;\"> 0.676 </td>\n   <td style=\"text-align:center;\"> 0.933 </td>\n   <td style=\"text-align:center;\"> 0.371 </td>\n   <td style=\"text-align:center;\"> 0.042 </td>\n   <td style=\"text-align:center;\"> 0.423 </td>\n   <td style=\"text-align:center;\"> 0.716 </td>\n   <td style=\"text-align:center;\"> 0.147 </td>\n  </tr>\n  <tr>\n   \n   <td style=\"text-align:center;vertical-align: middle !important;\" rowspan=\"2\"> ZM96.gp140 </td>\n   <td style=\"text-align:center;\"> 6.5 </td>\n   <td style=\"text-align:center;\"> 0.030 </td>\n   <td style=\"text-align:center;\"> 0.631 </td>\n   <td style=\"text-align:center;\"> 0.244 </td>\n   <td style=\"text-align:center;\"> 0.750 </td>\n   <td style=\"text-align:center;\"> 0.234 </td>\n   <td style=\"text-align:center;\"> 0.949 </td>\n   <td style=\"text-align:center;\"> 0.890 </td>\n   <td style=\"text-align:center;\"> 0.994 </td>\n   <td style=\"text-align:center;\"> 0.090 </td>\n   <td style=\"text-align:center;\"> 0.504 </td>\n  </tr>\n  <tr>\n   \n   \n   <td style=\"text-align:center;\"> 12.0 </td>\n   <td style=\"text-align:center;\"> 0.186 </td>\n   <td style=\"text-align:center;\"> 0.353 </td>\n   <td style=\"text-align:center;\"> 0.286 </td>\n   <td style=\"text-align:center;\"> 0.722 </td>\n   <td style=\"text-align:center;\"> 0.889 </td>\n   <td style=\"text-align:center;\"> 0.409 </td>\n   <td style=\"text-align:center;\"> 0.459 </td>\n   <td style=\"text-align:center;\"> 0.376 </td>\n   <td style=\"text-align:center;\"> 0.021 </td>\n   <td style=\"text-align:center;\"> 0.189 </td>\n  </tr>\n  <tr>\n   <td style=\"text-align:center;vertical-align: middle !important;\" rowspan=\"2\"> CD4+ </td>\n   <td style=\"text-align:center;vertical-align: middle !important;\" rowspan=\"2\"> Any ENV PTEG </td>\n   <td style=\"text-align:center;\"> 6.5 </td>\n   <td style=\"text-align:center;\"> 0.237 </td>\n   <td style=\"text-align:center;\"> 0.555 </td>\n   <td style=\"text-align:center;\"> 0.771 </td>\n   <td style=\"text-align:center;\"> 0.108 </td>\n   <td style=\"text-align:center;\"> 0.237 </td>\n   <td style=\"text-align:center;\"> 0.104 </td>\n   <td style=\"text-align:center;\"> 0.717 </td>\n   <td style=\"text-align:center;\"> 0.084 </td>\n   <td style=\"text-align:center;\"> 0.627 </td>\n   <td style=\"text-align:center;\"> 0.319 </td>\n  </tr>\n  <tr>\n   \n   \n   <td style=\"text-align:center;\"> 12.0 </td>\n   <td style=\"text-align:center;\"> 0.248 </td>\n   <td style=\"text-align:center;\"> 0.565 </td>\n   <td style=\"text-align:center;\"> 0.434 </td>\n   <td style=\"text-align:center;\"> 0.103 </td>\n   <td style=\"text-align:center;\"> 0.117 </td>\n   <td style=\"text-align:center;\"> 0.856 </td>\n   <td style=\"text-align:center;\"> 0.861 </td>\n   <td style=\"text-align:center;\"> 0.128 </td>\n   <td style=\"text-align:center;\"> 0.773 </td>\n   <td style=\"text-align:center;\"> 0.723 </td>\n  </tr>\n</tbody>\n</table>"

allMDS.html %>% cat(file = 'tables/allMDS.html')

allMDS %>%
kable("latex", escape = F, digits = 3, align = 'c', booktabs = T) %>%
collapse_rows(columns = 1:2, latex_hline = "full") %>%
pass -> allMDS.latex


allMDS.latex %>% cat(file = 'tables/allMDS.tex')

# Print latex table to tex file
cat(docHead, allMDS.latex, docTail, file = 'tables/allMDS.tex')

write_csv(allMDS, 'tables/allMDSGlmPValues.csv')

Agglomeration

# How many taxa do we see if we agglomerate at different levels
psN2 %>% tax_table %>% as.data.frame %>% dplyr::select(Phylum:Genus) %>% colnames -> taxLevels

data_frame(taxLevels) %>%
mutate(ntaxa = map(taxLevels,
    function(lev){
        psN2 %>% tax_glom(lev) %>% ntaxa
    }
                                             )) %>%
mutate(ntaxa = unlist(ntaxa)) %>%
pass -> NTaxaAtLevel
NTaxaAtLevel

## # A tibble: 5 x 2
##   taxLevels ntaxa
##   <chr>     <int>
## 1 Phylum        5
## 2 Class        12
## 3 Order        17
## 4 Family       39
## 5 Genus       104

data_frame(taxLevels = "Species", ntaxa = ntaxa(psN2), ps = list((psN2))) -> specRow
data_frame(taxLevels = "Species", ntaxa = ntaxa(psN1), psCount = list((psN1))) -> specRowC

D2K_savename <- function(distmat){
    # cascade names forward with the D2K operation
    require(MiRKAT)
    out <- MiRKAT::D2K(distmat)
    colnames(out) <- colnames(distmat)
    rownames(out) <- rownames(distmat)
    out
}

# Data frame of phyloseq objects distances and kernels at a bunch of taxonomic levels
NTaxaAtLevel %>%
mutate(ps = map(ntaxa, ~tip_glom_saveid(psN2, k = .))) %>%
# process the phyloseq objects so they have better names
mutate(ps = map(ps, ~swap.phyloseq.taxnames(tag_phyloseq(remove_tag_phyloseq(.)), oldname = 'oldname2'))) %>%
# add in the species data row (which should already have correct names)
bind_rows(specRow) %>%
# calculate jensen-shannon distance matrix
mutate(jsd = map(ps, ~phyloseq::distance(., method = "jsd") )) %>%
# convert to 2d matrix
mutate(jsdMat = map(jsd, ~as.matrix(.))) %>%
# calculate kernel
mutate(kjsd = map(jsdMat, ~D2K_savename(.))) -> tmp

## Loading required package: cluster

## Warning: Setting row names on a tibble is deprecated.

## Warning: Setting row names on a tibble is deprecated.

## Warning: Setting row names on a tibble is deprecated.

## Warning: Setting row names on a tibble is deprecated.

## Warning: Setting row names on a tibble is deprecated.

tmp %>%
mutate(psNoZero = map(ps, ~transform_sample_counts(., function(x) x+(1/1000)))) -> tmp

tmp %>%
## chemometrics::clr just works, while compositions::clr throws a criptic error message here
mutate(clr = map(psNoZero, ~ transform_otu_table(., chemometrics::clr))) %>%
#mutate(clr = map(psNoZero, ~ transform_otu_table(., function(x) as.matrix(compositions::clr(x))))) %>%
pass -> psDf0 # Original way

# Data frame of phyloseq objects distances and kernels at a bunch of taxonomic levels
# I use psN1 because I need count data for some downstream steps.
NTaxaAtLevel %>%
mutate(psCount = map(ntaxa, ~tip_glom_saveid(psN1, k = .))) %>%
# process the phyloseq objects so they have better names
mutate(psCount = map(psCount, ~swap.phyloseq.taxnames(tag_phyloseq(remove_tag_phyloseq(.)), oldname = 'oldname2'))) %>%
# add in the species data row (which should already have correct names)
bind_rows(specRowC) %>%
pass -> tmp

## Warning: Setting row names on a tibble is deprecated.

## Warning: Setting row names on a tibble is deprecated.

## Warning: Setting row names on a tibble is deprecated.

## Warning: Setting row names on a tibble is deprecated.

## Warning: Setting row names on a tibble is deprecated.

tmp %>%
# calculate jensen-shannon distance matrix
mutate(ps = map(psCount, ~transform_sample_counts(., function(x) {x/sum(x)}))) %>%
mutate(jsd = map(ps, ~phyloseq::distance(., method = "jsd") )) %>%
# convert to 2d matrix
mutate(jsdMat = map(jsd, ~as.matrix(.))) %>%
# calculate kernel
mutate(kjsd = map(jsdMat, ~D2K_savename(.))) -> tmp

tmp %>%
mutate(psNoZero = map(ps, ~transform_sample_counts(., function(x) x+(1/1000)))) %>%
## chemometrics::clr just works, while compositions::clr throws a criptic error message here
mutate(clr = map(psNoZero, ~ transform_otu_table(., chemometrics::clr))) %>%
#mutate(clr = map(psNoZero, ~ transform_otu_table(., function(x) as.matrix(compositions::clr(x))))) %>%
pass -> psDf

print(psDf)

## # A tibble: 6 x 9
##   taxLevels ntaxa psCount   ps      jsd    jsdMat   kjsd   psNoZero clr   
##   <chr>     <int> <list>    <list>  <list> <list>   <list> <list>   <list>
## 1 Phylum        5 <S4: phy… <S4: p… <S3: … <dbl [2… <dbl … <S4: ph… <S4: …
## 2 Class        12 <S4: phy… <S4: p… <S3: … <dbl [2… <dbl … <S4: ph… <S4: …
## 3 Order        17 <S4: phy… <S4: p… <S3: … <dbl [2… <dbl … <S4: ph… <S4: …
## 4 Family       39 <S4: phy… <S4: p… <S3: … <dbl [2… <dbl … <S4: ph… <S4: …
## 5 Genus       104 <S4: phy… <S4: p… <S3: … <dbl [2… <dbl … <S4: ph… <S4: …
## 6 Species     651 <S4: phy… <S4: p… <S3: … <dbl [2… <dbl … <S4: ph… <S4: …

MirMulti <- function(x, KsDf = psDf, ps = psN2, nperm = 9999){
    
    Ks = KsDf$kjsd
    
    # I  bind to the phyloseq object and then peel off again later to guerentee
    # that the y-data is in the same order as the Ks
    locPS <- phylo_join(ps, x, by = 'pub_id')
    
    ydata0 <- sample_data(locPS)$mag
    yna <- is.na(ydata0)

    ydata <- ydata0[!yna]
    loc.Ks <- lapply(Ks, function(K){K[!yna, !yna]})  
  
    bcxJSD <- MiRKAT(y = jac_box_cox_2(ydata), Ks = loc.Ks, out_type = "C", method = 'permutation', nperm = nperm)
    medJSD <- MiRKAT(y = medcode(ydata), Ks = loc.Ks, out_type = "D", method = 'permutation', nperm = nperm)
    mmDf = data.frame(
        taxLevels = KsDf$taxLevels,
        ntaxa = KsDf$ntaxa,
        bcxJSD = bcxJSD$indivP, medJSD = medJSD$indivP,
        bcxJSDOmni = bcxJSD$omnibus_p, medJSDOmni = medJSD$omnibus_p)
    mmDf
    
    }

# Test case

use.immune %>%
filter(type == 'IgG' & antigen == 'gp41'& visitno == 2 & ct == 'T') -> test.immune1

test.mm <- MirMulti(test.immune1, Ks = psDf, nperm = 999)

test.mm

##   taxLevels ntaxa bcxJSD medJSD bcxJSDOmni medJSDOmni
## 1    Phylum     5  0.217  0.228      0.162     0.0435
## 2     Class    12  0.282  0.228      0.162     0.0435
## 3     Order    17  0.048  0.011      0.162     0.0435
## 4    Family    39  0.123  0.031      0.162     0.0435
## 5     Genus   104  0.247  0.038      0.162     0.0435
## 6   Species   651  0.505  0.131      0.162     0.0435

ptm = proc.time()

use.immune %>%
group_by(type, antigen, month) %>%
nest %>%
mutate(mir = map(data,
    ~MirMulti(., Ks = psDf, ps = psN2, nperm = 999)
)) %>%
dplyr::select(-data) %>% unnest(mir) %>%
pass -> mirLevels

proc.time() - ptm

##    user  system elapsed 
##   2.561   0.000   2.561

mirLevels %>% dplyr::select(-ntaxa, -medJSD) %>% spread(key = taxLevels, value = bcxJSD)

## # A tibble: 20 x 11
##    type  antigen month bcxJSDOmni medJSDOmni  Class Family  Genus  Order
##    <chr> <chr>   <dbl>      <dbl>      <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
##  1 CD4+  ANY.EN…   6.5     0.762      0.285  0.497  0.512  0.704  0.396 
##  2 CD4+  ANY.EN…  12       0.054      0.148  0.014  0.146  0.136  0.057 
##  3 IgA   gp41      0       0.995      0.864  0.997  0.875  0.916  0.953 
##  4 IgA   gp41      6.5     0.259      0.362  0.316  0.546  0.401  0.395 
##  5 IgA   gp41     12       0.525      0.916  0.878  0.419  0.349  0.712 
##  6 IgA   p24       0       0.886      0.388  0.921  0.815  0.897  0.917 
##  7 IgA   p24       6.5     0.982      0.830  0.905  0.933  0.8    0.852 
##  8 IgA   p24      12       0.444      0.102  0.167  0.604  0.639  0.298 
##  9 IgG   Con.6.…   6.5     0.388      0.052  0.233  0.178  0.178  0.208 
## 10 IgG   Con.6.…  12       0.003      0.053  0.04   0.002  0.001  0.002 
## 11 IgG   gp41      0       0.158      0.0365 0.269  0.124  0.229  0.047 
## 12 IgG   gp41      6.5     0.334      0.0765 0.115  0.26   0.148  0.131 
## 13 IgG   gp41     12       0.310      0.722  0.156  0.3    0.179  0.135 
## 14 IgG   gp70_B…   6.5     0.526      0.949  0.217  0.415  0.429  0.376 
## 15 IgG   gp70_B…  12       0.114      0.086  0.169  0.0720 0.0920 0.032 
## 16 IgG   p24       0       0.893      0.476  0.718  0.951  0.611  0.913 
## 17 IgG   p24       6.5     0.0185     0.853  0.005  0.051  0.028  0.047 
## 18 IgG   p24      12       0.782      0.696  0.971  0.469  0.558  0.41  
## 19 IgG   ZM96.g…   6.5     0.097      0.0405 0.0860 0.113  0.057  0.124 
## 20 IgG   ZM96.g…  12       0.280      0.350  0.464  0.184  0.105  0.0960
## # ... with 2 more variables: Phylum <dbl>, Species <dbl>

mirLevels %>% dplyr::select(-ntaxa, -bcxJSD) %>% spread(key = taxLevels, value = medJSD)

## # A tibble: 20 x 11
##    type  antigen month bcxJSDOmni medJSDOmni  Class Family  Genus  Order
##    <chr> <chr>   <dbl>      <dbl>      <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
##  1 CD4+  ANY.EN…   6.5     0.762      0.285  0.102   0.193 0.268  0.0970
##  2 CD4+  ANY.EN…  12       0.054      0.148  0.047   0.164 0.184  0.061 
##  3 IgA   gp41      0       0.995      0.864  0.855   0.824 0.905  0.917 
##  4 IgA   gp41      6.5     0.259      0.362  0.652   0.485 0.43   0.364 
##  5 IgA   gp41     12       0.525      0.916  0.918   0.683 0.647  0.818 
##  6 IgA   p24       0       0.886      0.388  0.92    0.559 0.791  0.847 
##  7 IgA   p24       6.5     0.982      0.830  0.845   0.831 0.694  0.556 
##  8 IgA   p24      12       0.444      0.102  0.028   0.24  0.401  0.114 
##  9 IgG   Con.6.…   6.5     0.388      0.052  0.0830  0.019 0.014  0.015 
## 10 IgG   Con.6.…  12       0.003      0.053  0.0990  0.038 0.042  0.015 
## 11 IgG   gp41      0       0.158      0.0365 0.211   0.03  0.032  0.01  
## 12 IgG   gp41      6.5     0.334      0.0765 0.059   0.06  0.0690 0.021 
## 13 IgG   gp41     12       0.310      0.722  0.905   0.747 0.332  0.698 
## 14 IgG   gp70_B…   6.5     0.526      0.949  0.723   0.949 0.974  0.906 
## 15 IgG   gp70_B…  12       0.114      0.086  0.525   0.044 0.025  0.023 
## 16 IgG   p24       0       0.893      0.476  0.397   0.287 0.182  0.752 
## 17 IgG   p24       6.5     0.0185     0.853  0.479   0.673 0.502  0.69  
## 18 IgG   p24      12       0.782      0.696  0.745   0.329 0.374  0.348 
## 19 IgG   ZM96.g…   6.5     0.097      0.0405 0.0820  0.062 0.0720 0.043 
## 20 IgG   ZM96.g…  12       0.280      0.350  0.143   0.523 0.498  0.341 
## # ... with 2 more variables: Phylum <dbl>, Species <dbl>

I’d like to combine the above into one table, when it isn’t 7:45. Probably has soemething to do with merging columns or something. Or maybe I just want to plot it as a figure.

mirLevels %>%
gather(metric, P, bcxJSD:medJSD) -> mirDat

mirDat %>% dplyr::select(type:month, bcxJSD = bcxJSDOmni, medJSD = medJSDOmni) %>%
group_by(type, antigen, month) %>%
summarize(bcxJSD = mean(bcxJSD), medJSD = mean(medJSD)) %>%
gather(metric, P, bcxJSD, medJSD) -> mirOmni

mirOmni

## # A tibble: 40 x 5
## # Groups:   type, antigen [8]
##    type  antigen       month metric     P
##    <chr> <chr>         <dbl> <chr>  <dbl>
##  1 CD4+  ANY.ENV.PTEG    6.5 bcxJSD 0.762
##  2 CD4+  ANY.ENV.PTEG   12   bcxJSD 0.054
##  3 IgA   gp41            0   bcxJSD 0.995
##  4 IgA   gp41            6.5 bcxJSD 0.259
##  5 IgA   gp41           12   bcxJSD 0.525
##  6 IgA   p24             0   bcxJSD 0.886
##  7 IgA   p24             6.5 bcxJSD 0.982
##  8 IgA   p24            12   bcxJSD 0.444
##  9 IgG   Con.6.gp120.B   6.5 bcxJSD 0.388
## 10 IgG   Con.6.gp120.B  12   bcxJSD 0.003
## # ... with 30 more rows

NTaxaAtLevel %>% bind_rows(specRow[,1:2]) %>% unite(nLev, taxLevels, ntaxa, remove = FALSE) -> NTaxaAtLevel2
NTaxaAtLevel2

## # A tibble: 6 x 3
##   nLev        taxLevels ntaxa
##   <chr>       <chr>     <int>
## 1 Phylum_5    Phylum        5
## 2 Class_12    Class        12
## 3 Order_17    Order        17
## 4 Family_39   Family       39
## 5 Genus_104   Genus       104
## 6 Species_651 Species     651

bind_rows(
    permKernTable %>% mutate(metric = 'med'),
    permKernTableGaus %>% mutate(metric = 'bcx')
    ) %>%
dplyr::select(type, antigen, month, metric, mir.P) %>%
pass -> WufPData

fixant <- function(df){
    df %>%
    mutate(antigen = stringr::str_replace_all(antigen, "\\.", " ")) %>%
    mutate(antigen = stringr::str_replace_all(antigen, "_", " ")) %>%
    #mutate(metric =  stringr::str_replace_all(metric, "bcx", "")) %>%
    pass
}

fixstuff <- function(df){
    df %>%
    fixant %>%
    mutate(metric =  stringr::str_replace_all(metric, "JSD", "")) %>%
    pass
}

mirDat %>% fixant %>% head

## # A tibble: 6 x 9
##   type  antigen month taxLevels ntaxa bcxJSDOmni medJSDOmni metric      P
##   <chr> <chr>   <dbl> <fct>     <int>      <dbl>      <dbl> <chr>   <dbl>
## 1 IgA   gp41      6.5 Phylum        5      0.259      0.362 bcxJSD 0.0970
## 2 IgA   gp41      6.5 Class        12      0.259      0.362 bcxJSD 0.316 
## 3 IgA   gp41      6.5 Order        17      0.259      0.362 bcxJSD 0.395 
## 4 IgA   gp41      6.5 Family       39      0.259      0.362 bcxJSD 0.546 
## 5 IgA   gp41      6.5 Genus       104      0.259      0.362 bcxJSD 0.401 
## 6 IgA   gp41      6.5 Species     651      0.259      0.362 bcxJSD 0.731

mirDat %>% 
#mutate(antigen = stringr::str_replace_all(antigen, "\\.", " ")) %>% 
fixstuff %>%
ggplot(aes(x = ntaxa, y = P, col = factor(month), fill = factor(month))) +
geom_point(pch = 21) +
facet_grid(type + antigen ~ metric, labeller = labeller(antigen = label_wrap_gen(width = 10))) + 
scale_x_log10(breaks = c(3, NTaxaAtLevel2$ntaxa, 1000), labels = c("omni", NTaxaAtLevel2$nLev, "wunifrac")) +
scale_y_log10(breaks = c(0.002, 0.01, 0.05, 0.2, 1)) + 
geom_hline(yintercept=0.05, col = 'blue', alpha = 0.5) + geom_hline(yintercept=0.01, col = 'red', alpha = 0.5) +
#geom_hline(data = mirOmni, aes(yintercept = P, col = factor(month))) +
#annotation_logticks(sides = 'bl') +
#geom_rug(data = mirOmni, aes(y = P, col = factor(month)), inherit.aes = F) +
geom_point(data = mirOmni %>% ungroup %>% fixstuff,
           aes(x = 3, y = P, col = factor(month), fill = factor(month)), inherit.aes = F, pch = 22, size = 2) +
geom_point(data = WufPData %>% ungroup %>% fixant,
           aes(x = 1000, y = mir.P, col = factor(month), fill = factor(month)), inherit.aes = F, pch = 24, size = 2) +

scale_colour_manual(values=cbPalette) + 
scale_fill_manual(values=alpha(cbPalette, 0.5)) + 
labs(col = "month", fill = "month") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, hjust = 1,, vjust = 0.5, size = 10),
     strip.text.y = (element_text(angle = 90))) -> pjsd0

options(repr.plot.width=8, repr.plot.height= 6)
pjsd0

options(par0)
# I'd like to add weighted unifrac as a tick mark on the right.

NTaxaAtLevel2

## # A tibble: 6 x 3
##   nLev        taxLevels ntaxa
##   <chr>       <chr>     <int>
## 1 Phylum_5    Phylum        5
## 2 Class_12    Class        12
## 3 Order_17    Order        17
## 4 Family_39   Family       39
## 5 Genus_104   Genus       104
## 6 Species_651 Species     651

# New combined data frame that has omnibus, regular, and wunifrac all in one
bind_rows(
     mirDat %>% mutate(metric =  stringr::str_replace_all(metric, "JSD", "")) %>%
    mutate(test = "JSD"),
     mirOmni %>% ungroup %>% mutate(metric =  stringr::str_replace_all(metric, "JSD", "")) %>%
    mutate(taxLevels = "Omnibus") %>% mutate(test = "Omnibus"),
     WufPData %>% ungroup %>% dplyr::rename(P = mir.P) %>%
    mutate(taxLevels = "WUnifrac") %>% mutate(test = "WUnifrac")
) %>% 
mutate(antigen = factor(antigen, levels = c(ants2, ants1, "ANY.ENV.PTEG"))) %>%
mutate(type = factor(type, levels = c("IgA", "IgG", "CD4+"))) %>%
mutate(taxLevels = factor(taxLevels, levels = c("Omnibus", NTaxaAtLevel2$taxLevels, "WUnifrac"))) %>%
dplyr::select(-c(bcxJSDOmni:medJSDOmni))%>%
unite(nLev, taxLevels, ntaxa, remove = FALSE) %>%
mutate(nLev = stringr::str_replace(nLev, "_NA", "")) %>%
#mutate(antigen = stringr::str_replace_all(antigen, "\\.", " ")) %>%
#mutate(antigen = stringr::str_replace_all(antigen, "_", " ")) %>%
 mutate(antigen = factor(antigen, labels = stringr::str_replace_all(levels(antigen), "\\.", " "))) %>%
 mutate(antigen = factor(antigen, labels = stringr::str_replace_all(levels(antigen), "_", " "))) %>%


# mutate(antigen = factor(antigen, labels = (levels(antigen)))) %>%
# mutate(antigen = factor(antigen, labels = (levels(antigen)))) %>%

mutate(test = factor(test, levels = c('Omnibus', 'JSD', 'WUnifrac'))) %>%
pass -> mirDat2

## Warning in bind_rows_(x, .id): binding factor and character vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

mirDat2 %>% filter(type == 'CD4+')

## # A tibble: 32 x 9
##    type  antigen      month nLev        taxLevels ntaxa metric     P test 
##    <fct> <fct>        <dbl> <chr>       <fct>     <int> <chr>  <dbl> <fct>
##  1 CD4+  ANY ENV PTEG   6.5 Phylum_5    Phylum        5 bcx    0.86  JSD  
##  2 CD4+  ANY ENV PTEG   6.5 Class_12    Class        12 bcx    0.497 JSD  
##  3 CD4+  ANY ENV PTEG   6.5 Order_17    Order        17 bcx    0.396 JSD  
##  4 CD4+  ANY ENV PTEG   6.5 Family_39   Family       39 bcx    0.512 JSD  
##  5 CD4+  ANY ENV PTEG   6.5 Genus_104   Genus       104 bcx    0.704 JSD  
##  6 CD4+  ANY ENV PTEG   6.5 Species_651 Species     651 bcx    0.919 JSD  
##  7 CD4+  ANY ENV PTEG  12   Phylum_5    Phylum        5 bcx    0.908 JSD  
##  8 CD4+  ANY ENV PTEG  12   Class_12    Class        12 bcx    0.014 JSD  
##  9 CD4+  ANY ENV PTEG  12   Order_17    Order        17 bcx    0.057 JSD  
## 10 CD4+  ANY ENV PTEG  12   Family_39   Family       39 bcx    0.146 JSD  
## # ... with 22 more rows

mirDat2 %>% head

## # A tibble: 6 x 9
##   type  antigen month nLev        taxLevels ntaxa metric      P test 
##   <fct> <fct>   <dbl> <chr>       <fct>     <int> <chr>   <dbl> <fct>
## 1 IgA   gp41      6.5 Phylum_5    Phylum        5 bcx    0.0970 JSD  
## 2 IgA   gp41      6.5 Class_12    Class        12 bcx    0.316  JSD  
## 3 IgA   gp41      6.5 Order_17    Order        17 bcx    0.395  JSD  
## 4 IgA   gp41      6.5 Family_39   Family       39 bcx    0.546  JSD  
## 5 IgA   gp41      6.5 Genus_104   Genus       104 bcx    0.401  JSD  
## 6 IgA   gp41      6.5 Species_651 Species     651 bcx    0.731  JSD

mirDat2 %>%
filter(type == "IgG") %>%
ggplot(aes(x = taxLevels, y = P, col = factor(month), fill = factor(month), shape = test)) +
geom_point(size = 2) +
facet_grid(type + antigen ~ metric, labeller = labeller(antigen = label_wrap_gen(width = 10))) + 
#scale_x_log10(breaks = c(3, NTaxaAtLevel2$ntaxa, 1000), labels = c("omni", NTaxaAtLevel2$nLev, "wunifrac")) +
scale_y_log10(breaks = c(0.002, 0.01, 0.05, 0.2, 1)) + 
geom_hline(yintercept=0.05, col = 'blue', alpha = 0.5) + geom_hline(yintercept=0.01, col = 'red', alpha = 0.5) +

scale_shape_manual(values = c(22, 21, 24)) +
scale_colour_manual(values=cbPalette) + 
scale_fill_manual(values=alpha(cbPalette, 0.5)) + 
labs(col = "month", fill = "month", x = "Taxonomic Level") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, hjust = 1,, vjust = 0.5, size = 10),
     strip.text.y = (element_text(angle = 90))) -> pjsd


options(repr.plot.width=6, repr.plot.height= 8)
pjsd

ggsave('figures/KernelPVsLevel.png', width = 6, height = 8)
options(par0)

X-axis is now spaced evenly

Table SX. P values of kernel regression tests. Circles indicate jensen shannon values at different taxonomic resolutions. Squares are the omnibus p-value for that cohort of tests. Triangles indicate kernel regression p-values for the corresponding weighted unifrac test.

The blue and red lines indicate p values of 5% and 1% respectively.

Observations: The weighted unifrac test is sensitive. In cases where only one taxonic level hits, weighted unifrac often also falls at some statistically significant value. The omnibus p value is often higher than the weighted unifrac one. Weighted unifrac seems like a good test for identifying patterns at any level that relate to an outcome. The jensen shannon informs us about which level the pattern is observed.

Family, genera and species vs wuf1

I might even be able to drill down to every level.

model_each_species <- function(ps, f, pthresh = 1, q = FALSE){
    # Start with the otu table
ps %>%
# reshape it so we have clr values for every taxon-sample pair
otu_table %>% as.data.frame %>% rownames_to_column("Sample") %>% gather(Taxon, clr, -Sample) %>%
    # bind that to the sample data
    # doing this here seems remarkably inefficient, but its not creating a bottleneck so I'll leave it.
left_join(
    ps %>%
    # the sample data need to have MDS1 and MDS2 appended to them
    phylo_join(
    psN2.pcoa %>% scores(display = "sites") %>% # hardcoded psN2.pcoa
        as.data.frame %>% 
        rownames_to_column %>% 
        dplyr::select('rowname', 'MDS1', 'MDS2'),
    by = 'rowname'
) %>%
    # back to binding to sample data
    sample_data %>% as('data.frame') %>% rownames_to_column("Sample"),
     by = 'Sample') %>%

group_by(Taxon) %>%  # group and nest for model run
nest %>%
mutate(Mod = map(data, f)) %>% # apply model over each species
mutate(Glance = map(Mod, glance), Tidy = map(Mod, tidy)) %>% # extract relevant data from model
# view model
dplyr::select(Taxon, Tidy) %>% unnest %>%
mutate(term = gsub('[\\( \\)]','', term)) %>% # remove parentheses from "(Intercept)"
gather(meas, val, estimate:p.value) %>% 
unite(meas, term, meas) %>% spread(meas, val) %>% arrange(clr_estimate) %>% 
dplyr::select(Taxon, Intercept_estimate, clr_estimate, clr_std.error, clr_p.value) %>%
    # add q value
    {if(q) mutate(., clr_q.value = p2q(clr_p.value)) else .} %>%
    
 filter(clr_p.value < pthresh) %>%

     #Join taxonomy information
     left_join(
     as.data.frame(ps@tax_table@.Data) %>% as.tibble %>% dplyr::select(Kingdom:Genus, Species, tag) %>%
         mutate(tag = as.character(tag)), # mutate so tag is and taxon are both character class
     by = c("Taxon" = "tag")) %>%
pass
 }

model_each_species_for_antigen <- function(antigen, ps = psN2){
    ps %>%
    model_each_species(function(df){glm(medcode2(get(antigen)) ~ clr, data = df, family = 'binomial')}, q = TRUE, pthresh = 1)
}

ColsToRun <- c('IgG_Con.6.gp120.B_Month_6.5', 'IgG_Con.6.gp120.B_Month_12', 'IgG_gp41_Month_0', 'IgG_gp41_Month_6.5', 'IgG_gp70_B.CaseA_V1_V2_Month_12', 'IgG_ZM96.gp140_Month_6.5', 'MDS1', 'isMale' ) 

model_each_species_case <- function(ps){
    
    ps %>% model_each_species(function(df){glm(MDS1 ~ clr, data = df, family = 'gaussian')}, q = TRUE, pthresh = 1) %>%
arrange(clr_estimate) %>%
mutate(Taxon = factor(Taxon, levels = Taxon[order(clr_estimate)])) %>%
    mutate(test = 'gaussian', antigen = 'MDS1') %>%
pass -> loc_mds1Glms
    
        ps %>% model_each_species(function(df){glm(log10(IgG_Con.6.gp120.B_Month_12 + 100) ~ clr, data = df, family = 'gaussian')}, q = TRUE, pthresh = 1) %>%
arrange(clr_estimate) %>%
 mutate(Taxon = factor(Taxon, levels = levels(loc_mds1Glms$Taxon))) %>%
    mutate(test = 'gaussian', antigen = 'Con.6.gp120.B_Month_12') %>%
 pass -> loc_gp120Glms
    
      tibble(antigen = ColsToRun) %>% mutate(model = map(antigen, ~model_each_species_for_antigen(., ps = ps))) %>%
  unnest %>% mutate(Taxon = factor(Taxon, levels = levels(loc_mds1Glms$Taxon))) %>%
        mutate(test = 'binomial') %>%
    pass-> loc_logitCoefs
    

    #list(loc_mds1Glms, loc_gp120Glms, loc_logitCoefs)
     bind_rows(loc_mds1Glms, loc_gp120Glms, loc_logitCoefs) %>% dplyr::select(test, antigen, everything())
    
}

#psDf %>% mutate(ps2 = map(ps, ~swap.phyloseq.taxnames(tag_phyloseq(remove_tag_phyloseq(.))))) -> test

#psDf[[1,"clr"]] %>% tax_table

psDf[[1]]

## [1] "Phylum"  "Class"   "Order"   "Family"  "Genus"   "Species"

psDf %>% print

## # A tibble: 6 x 9
##   taxLevels ntaxa psCount   ps      jsd    jsdMat   kjsd   psNoZero clr   
##   <chr>     <int> <list>    <list>  <list> <list>   <list> <list>   <list>
## 1 Phylum        5 <S4: phy… <S4: p… <S3: … <dbl [2… <dbl … <S4: ph… <S4: …
## 2 Class        12 <S4: phy… <S4: p… <S3: … <dbl [2… <dbl … <S4: ph… <S4: …
## 3 Order        17 <S4: phy… <S4: p… <S3: … <dbl [2… <dbl … <S4: ph… <S4: …
## 4 Family       39 <S4: phy… <S4: p… <S3: … <dbl [2… <dbl … <S4: ph… <S4: …
## 5 Genus       104 <S4: phy… <S4: p… <S3: … <dbl [2… <dbl … <S4: ph… <S4: …
## 6 Species     651 <S4: phy… <S4: p… <S3: … <dbl [2… <dbl … <S4: ph… <S4: …

ptm = proc.time()
psDf %>% dplyr::select(taxLevels, ntaxa, clr) %>% mutate(localmod = map(clr, model_each_species_case)) ->psDfLoc

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

proc.time() - ptm

##    user  system elapsed 
##  41.692   0.032  41.727

print(psDfLoc)

## # A tibble: 6 x 4
##   taxLevels ntaxa clr            localmod             
##   <chr>     <int> <list>         <list>               
## 1 Phylum        5 <S4: phyloseq> <tibble [50 × 15]>   
## 2 Class        12 <S4: phyloseq> <tibble [120 × 15]>  
## 3 Order        17 <S4: phyloseq> <tibble [170 × 15]>  
## 4 Family       39 <S4: phyloseq> <tibble [390 × 15]>  
## 5 Genus       104 <S4: phyloseq> <tibble [1,040 × 15]>
## 6 Species     651 <S4: phyloseq> <tibble [6,510 × 15]>

psDfLoc %>% dplyr::select(-clr) %>% unnest(localmod) -> tmp

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector

psDfLoc$taxLevels

## [1] "Phylum"  "Class"   "Order"   "Family"  "Genus"   "Species"

tmp %>% mutate(taxLevels = factor(taxLevels, levels = psDfLoc$taxLevels)) -> LocalTests

LocalTests %>% 
filter(antigen != "MDS1") %>%
write_csv("tables/AllLocalTests.csv")

I want to show the local tests vs antibodies.

LocalTests %>% head

## # A tibble: 6 x 17
##   taxLevels ntaxa test  antigen Taxon Intercept_estim… clr_estimate
##   <fct>     <int> <chr> <chr>   <chr>            <dbl>        <dbl>
## 1 Phylum        5 gaus… MDS1    Bact…            1.94        -0.762
## 2 Phylum        5 gaus… MDS1    Porp…           -1.03        -0.514
## 3 Phylum        5 gaus… MDS1    Acti…           -0.201       -0.104
## 4 Phylum        5 gaus… MDS1    Camp…            0.439        0.242
## 5 Phylum        5 gaus… MDS1    Bact…           -1.42         0.445
## 6 Phylum        5 gaus… Con.6.… Bact…            3.76        -0.285
## # ... with 10 more variables: clr_std.error <dbl>, clr_p.value <dbl>,
## #   clr_q.value <dbl>, Kingdom <fct>, Phylum <chr>, Class <chr>,
## #   Order <chr>, Family <chr>, Genus <chr>, Species <chr>

options(repr.plot.width=6, repr.plot.height=8)
LocalTests %>% 
filter(antigen != "MDS1") %>%
#Clean up labels
mutate(antigen = stringr::str_replace_all(antigen, "_", " ")) %>%
mutate(antigen = stringr::str_replace_all(antigen, " Month", " -- Month")) %>%
mutate(antigen = stringr::str_replace_all(antigen, " IgG ", "")) %>%
ggplot(aes(x = factor(taxLevels), y = clr_q.value)) + geom_point(size = 0.1) +
geom_hline(yintercept = 0.2, color = 'grey') +
 geom_hline(yintercept = 0.05, color = 'blue') + geom_hline(yintercept = 0.01, color = 'green') +
facet_wrap(~antigen + test) + scale_y_log10() + 
  labs(y = "q-value", x = "Taxonomic Level") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))

ggsave('figures/LocalQEveryLevel.png', width = 6, height = 8, units = "in")

BOOKMARK

options(repr.plot.width=6, repr.plot.height=8)
LocalTests %>% 
filter(antigen != "MDS1") %>%
#Clean up labels
mutate(antigen = stringr::str_replace_all(antigen, "_", " ")) %>%
mutate(antigen = stringr::str_replace_all(antigen, " Month", " -- Month")) %>%
mutate(antigen = stringr::str_replace_all(antigen, " IgG ", "")) %>%
ggplot(aes(x = factor(taxLevels), y = clr_p.value)) + geom_point(size = 0.1) +
geom_hline(yintercept = 0.2, color = 'grey') +
 geom_hline(yintercept = 0.05, color = 'blue') + geom_hline(yintercept = 0.01, color = 'green') +
facet_wrap(~antigen + test) +# scale_y_log10() +
  labs(y = expression(paste(italic("p"),"-value")), x = "Taxonomic Level") + 
  
theme_bw() +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))

ggsave('figures/LocalPEveryLevel.png', width = 6, height = 8, units = "in")

options(repr.plot.width=6, repr.plot.height=8)
LocalTests %>% 
filter(antigen != "MDS1") %>%
#Clean up labels
mutate(antigen = stringr::str_replace_all(antigen, "_", " ")) %>%
mutate(antigen = stringr::str_replace_all(antigen, " Month", " -- Month")) %>%
mutate(antigen = stringr::str_replace_all(antigen, " IgG ", "")) %>%
ggplot(aes(x = factor(taxLevels), y = clr_p.value)) + geom_point(size = 0.1) +
geom_hline(yintercept = 0.2, color = 'grey') +
 geom_hline(yintercept = 0.05, color = 'blue') + geom_hline(yintercept = 0.01, color = 'green') +
facet_wrap(~antigen + test) + scale_y_log10() +
labs(y = expression(paste(italic("p"),"-value")), x = "Taxonomic Level") + 
theme_bw() +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))

ggsave('figures/LocalPEveryLevel_LogScale.png', width = 6, height = 8, units = "in")

order_taxa_by_mds1 <- function(df){
    # this has to be a model_each_species type of data frame
    df %>% filter(antigen == 'MDS1' & test == 'gaussian') %>%
    mutate(TaxonF = factor(Taxon, levels = Taxon[order(clr_estimate)])) -> mds1df
    df %>% mutate(TaxonF = factor(Taxon, levels = levels(mds1df$TaxonF)))
}

To my annoyance, everything is labeled with IgG except gp120_12

LocalTests %>%
filter(taxLevels == 'Family') %>%
order_taxa_by_mds1 %>%

filter(
    (antigen %in% c('IgG_gp41_Month_0', 'IgG_gp41_Month_6.5', 'IgG_Con.6.gp120.B_Month_6.5'))|
    (antigen == 'Con.6.gp120.B_Month_12' & test == 'gaussian')
      ) %>%
mutate(antigen = factor(antigen,
                        levels = c('IgG_gp41_Month_0', 'IgG_gp41_Month_6.5',
                                   'IgG_Con.6.gp120.B_Month_6.5', 'Con.6.gp120.B_Month_12'))) %>%
filter(clr_p.value < 0.05 & clr_q.value < 0.2) %>%
dplyr::select(antigen:clr_estimate) %>%
dplyr::select(-Intercept_estimate) %>%
mutate(cordir = sign(clr_estimate)) %>%
pass

## # A tibble: 11 x 4
##    antigen                Taxon                clr_estimate cordir
##    <fct>                  <chr>                       <dbl>  <dbl>
##  1 Con.6.gp120.B_Month_12 Bacteroidetes.4            -0.307     -1
##  2 Con.6.gp120.B_Month_12 Firmicutes.4               -0.295     -1
##  3 Con.6.gp120.B_Month_12 Firmicutes.2               -0.157     -1
##  4 Con.6.gp120.B_Month_12 Clostridia                 -0.143     -1
##  5 Con.6.gp120.B_Month_12 Bacteroides                 0.204      1
##  6 Con.6.gp120.B_Month_12 Firmicutes.5                0.207      1
##  7 IgG_gp41_Month_0       Porphyromonadaceae.1       -2.12      -1
##  8 IgG_gp41_Month_0       Porphyromonadaceae         -1.69      -1
##  9 IgG_gp41_Month_0       Bacteroidetes.1            -1.20      -1
## 10 IgG_gp41_Month_0       Clostridia                  1.31       1
## 11 IgG_gp41_Month_0       Anaerococcus                1.53       1

LocalTests %>% pull(antigen) %>% unique

## [1] "MDS1"                            "Con.6.gp120.B_Month_12"         
## [3] "IgG_Con.6.gp120.B_Month_6.5"     "IgG_Con.6.gp120.B_Month_12"     
## [5] "IgG_gp41_Month_0"                "IgG_gp41_Month_6.5"             
## [7] "IgG_gp70_B.CaseA_V1_V2_Month_12" "IgG_ZM96.gp140_Month_6.5"       
## [9] "isMale"

labs(y = expression(paste(italic(“p”),“-value”)), x = “Taxonomic Level”) +

# Family Hits
LocalTests %>%
filter(taxLevels == 'Family') %>%
order_taxa_by_mds1 %>%

filter(
    (antigen %in% c('IgG_gp41_Month_0', 'IgG_gp41_Month_6.5', 'IgG_Con.6.gp120.B_Month_6.5', 
                    'IgG_ZM96.gp140_Month_6.5','IgG_gp70_B.CaseA_V1_V2_Month_12'))|
    (antigen == 'Con.6.gp120.B_Month_12' & test == 'gaussian')
      ) %>%
mutate(antigen = factor(antigen, levels = c(
    'IgG_gp41_Month_0', 'IgG_gp41_Month_6.5', 'IgG_Con.6.gp120.B_Month_6.5', 'Con.6.gp120.B_Month_12',
    'IgG_ZM96.gp140_Month_6.5','IgG_gp70_B.CaseA_V1_V2_Month_12'
))) %>%

pass -> tmp

#tmp$antigen %>% unique

tmp %>% filter(clr_p.value < 0.05 & clr_q.value < 0.2) %>%
pull(Taxon) %>% unique -> useFamily

tmp %>% filter(Taxon %in% useFamily) %>%

#Clean up labels
mutate(antigen = stringr::str_replace_all(antigen, "_", " ")) %>%
mutate(antigen = stringr::str_replace_all(antigen, " Month", " -- Month")) %>%
mutate(antigen = stringr::str_replace_all(antigen, "IgG ", "")) %>%

ggplot(aes(x = TaxonF, y = clr_estimate,
           color = (clr_p.value < 0.05), shape =(clr_q.value < 0.2))) +
geom_point(size = 3) + 
geom_errorbar(aes(ymin = clr_estimate - 2*clr_std.error, ymax = clr_estimate + 2*clr_std.error)) + 
  labs(y = "Coefficient",
       x = "Family Level Taxon",
       col = expression(paste(italic("p"), "-value")),
       shape = "q-value") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) + geom_hline(yintercept = 0) +
facet_wrap(~antigen, ncol = 1, scales = 'free_y')

# Show the censored ones accross - so this would be everything with at least one hit
# but also show what they are in all cases.

ggsave('figures/anyFamilyIgg.png', width = 6, height = 8)

I think its worth digging into clostridia and Prophyromonidaceae with stacked bars

Richness vs MDS1

btN_MDS <- betta(summary(bN1)$estimate,
              summary(bN1)$error,
              make_design_matrix(psN1, "MDS1")
)
btN_MDS

## $table
##             Estimates Standard Errors p-values
## (Intercept)  243.9720        16.82711    0.000
## predictors   -31.7198        27.31450    0.246
## 
## $cov
##             (Intercept) predictors
## (Intercept)    287.4724    57.2095
## predictors      57.2095   757.4669
## 
## $ssq_u
## [1] 1932.271
## 
## $homogeneity
## [1] 6.250059e+01 1.546382e-06
## 
## $global
## [1] 215.6907   0.0000
## 
## $blups
##  [1] 190.9549 235.6741 241.4954 196.3625 222.4246 270.5355 262.2978
##  [8] 280.2053 194.4011 282.1363 231.4248 227.0622 253.5892 260.3696
## [15] 266.7135 248.6283 296.0946 232.9051 201.2419 264.0743 264.8205

Not in any way that is statistically significant (p = 0.246)

Some pre-computing

richEsts <- bN1 %>% map_dbl("estimate")

richCI <- bN1 %>% map_df("interval") %>% t %>% as.data.frame() %>% rename(richLb = V1, richUb = V2) %>% merge(as.data.frame(richEsts), by = "row.names") %>% transform(row.names = Row.names, richEst = richEsts, Row.names = NULL, richEsts = NULL)
# need to add mds1 to this, or just tack this all onto psN1rich

psN1rich <- psN1
sample_data(psN1rich) <- merge(psN1@sam_data, richCI, by = "row.names") %>% transform(row.names = Row.names, Row.names = NULL)

Plot richness vs mds1

psN1rich %>% sample_data %>% as.data.frame %>%
  ggplot(aes(x = MDS1, y = richEst)) + geom_point() #+ geom_errorbar()

As above but with error bars.

psN1richP <- psN1rich %>% sample_data %>% as.data.frame %>%
  ggplot(aes(x = MDS1, y = richEst, ymin = richLb, ymax = richUb)) + geom_point() + geom_errorbar() + scale_y_log10()
psN1richP

Unifrac Distance vs Richness

The reviewer actually asked whether unifrac distance associates with alpha diversity. I’ve done that with MDS1, but not distance per se. Lets do one kernel regression, where we ask whether unifrac distance associates with richness.

Kernel regression, weighted unifrac vs richness

mirRich <- MiRKAT(y = richEsts, Ks = wufKN2, out_type = "C", method = 'permutation', nperm = jnperm)
mirRich

## [1] 0.22179

Very similar p-value to running betta against MDS1.

PCoA figure that compares MDS1 and MDS2 to richness

psN1richMDSP<- psN1rich %>%
sample_data() %>%
ggplot(aes(x = MDS1, y = MDS2)) + geom_point(aes(fill = richEst), size = 5, stroke = 1, shape = 21) +
viridis::scale_fill_viridis(name = 'richEst', direction = 1, option = "viridis") +
  coord_fixed(sqrt(psN2.pcoa$CA$eig[2]/psN2.pcoa$CA$eig[1])) +
#scale_colour_manual(name = 'gp41 Primary', values = c('black', 'grey70')) + 
cowplot::theme_cowplot()
psN1richMDSP

gpLocmtx <- cbind(1,
  psN1 %>% sample_data %>% as.matrix %>% as.data.frame %>% pull(IgG_gp41_Month_6.5) %>% as.character %>% as.numeric %>% medcode
)
colnames(gpLocmtx) = c("(Intercept)", "predictors")

btNLoc <- betta(summary(bN1)$estimate,
              summary(bN1)$error,
              gpLocmtx
)
btNLoc

## $table
##             Estimates Standard Errors p-values
## (Intercept) 209.10249        13.30004    0.000
## predictors   76.38459        22.30439    0.001
## 
## $cov
##             (Intercept) predictors
## (Intercept)    274.4924  -274.4924
## predictors    -274.4924   771.9782
## 
## $ssq_u
## [1] 780.9728
## 
## $homogeneity
## [1] 24.0326989  0.1949011
## 
## $global
## [1] 323.1194   0.0000
## 
## $blups
##  [1] 184.0879 256.1046 284.6409 195.8995 211.0920 290.4048 217.0720
##  [8] 288.8596 204.3579 291.4781 213.4344 210.8080 231.2472 222.1869
## [15] 291.5582 286.9768 296.3785 285.9490 200.8389 280.4231 287.5847

Compare richness to each immunogenicity parameter.

Step 1, make a data frame with BN1, but pub_id for each sample

Pre calculations

SampleNameToPubId <- psN1 %>% sample_data %>% as.data.frame %>% rownames_to_column() %>% as.tibble %>% dplyr::select(rowname, pub_id) %>% na.omit()

## Warning in class(x) <- c(subclass, "tbl_df", "tbl", "data.frame"): Setting
## class(x) to multiple strings ("tbl_df", "tbl", ...); result will no longer
## be an S4 object

A function that takes immunogenicity data, x (which we will pass in with tidyverse mapping functions), ad, the breakaway richness summary, and a transformation, defaults to medcode2 of the immunogenicity data.

immuneValpha <- function(x, ad = bN1, transformation = medcode2){
  #x <- arrange(x, rowname)
  x2 <- right_join(x, SampleNameToPubId, by = 'pub_id')
  x3 <- x2[is.finite(x2$mag),]
  
  ad2 <- ad[is.finite(x2$mag)]
  class(ad2) <- c("alpha_estimates", "list")
  
  gpXmtx = cbind(1, transformation(x3$mag))
  #gpXmtx
  
   thing <- betta(summary(ad2)$estimate,
               summary(ad2)$error,
               gpXmtx
 )
   out <- as.data.frame(t(thing$table[2,]))
   #colnames(out) = "pval"
   
   # add R^2 value, from simple pearson correlation
   locCor <- cor(summary(ad2)$estimate, gpXmtx[,2], method = "pearson")
   r2 <- locCor^2
   out2 <- data.frame(out, r2)
   
   out2
}

# test a single use case
immuneValpha(use.immune %>% filter(visitno == 9, type == "IgG", antigen == "Con.6.gp120.B"), bN1, transformation = medcode)

##   Estimates Standard.Errors p.values          r2
## 1  21.03933        22.41648    0.348 0.002631791

Actually run the analysis, raw table

immuneAlphaCompiled <- use.immune %>% 
  group_by(type, antigen, month) %>%
  do(immuneValpha(.)) %>%
  ungroup # if you forget to ungroup, pvalue calculations don't work correctly
immuneAlphaTable <- immuneAlphaCompiled %>% mutate(qval = p.adjust(`p.values`, method = "BH"))
immuneAlphaTable

## # A tibble: 20 x 8
##    type  antigen   month Estimates Standard.Errors p.values      r2   qval
##    <chr> <chr>     <dbl>     <dbl>           <dbl>    <dbl>   <dbl>  <dbl>
##  1 CD4+  ANY.ENV.…   6.5     40.7             19.9    0.041 3.26e-2 0.0911
##  2 CD4+  ANY.ENV.…  12       44.8             17.9    0.013 5.57e-3 0.0371
##  3 IgA   gp41        0       -1.84            27.3    0.946 3.22e-2 0.946 
##  4 IgA   gp41        6.5     16.3             26.2    0.534 1.45e-1 0.628 
##  5 IgA   gp41       12       68.4             28.7    0.017 2.11e-1 0.0425
##  6 IgA   p24         0        5.46            24.9    0.827 6.81e-3 0.871 
##  7 IgA   p24         6.5     56.6             21.9    0.01  1.60e-1 0.0333
##  8 IgA   p24        12       32.9             23.8    0.168 1.37e-2 0.28  
##  9 IgG   Con.6.gp…   6.5     21.0             22.4    0.348 2.63e-3 0.535 
## 10 IgG   Con.6.gp…  12       17.5             22.7    0.44  1.33e-3 0.587 
## 11 IgG   gp41        0       17.1             25.6    0.503 7.75e-4 0.628 
## 12 IgG   gp41        6.5     76.4             22.3    0.001 2.47e-1 0.01  
## 13 IgG   gp41       12       61.0             23.4    0.009 9.44e-2 0.0333
## 14 IgG   gp70_B.C…   6.5     -5.83            24.6    0.812 6.38e-3 0.871 
## 15 IgG   gp70_B.C…  12       70.5             14.7    0     6.21e-2 0     
## 16 IgG   p24         0      -30.8             20.3    0.129 3.33e-2 0.235 
## 17 IgG   p24         6.5     60.6             23.5    0.01  1.17e-1 0.0333
## 18 IgG   p24        12       47.3             26.4    0.074 2.34e-1 0.148 
## 19 IgG   ZM96.gp1…   6.5     21.9             25.3    0.387 3.55e-2 0.553 
## 20 IgG   ZM96.gp1…  12       66.8             22.0    0.002 9.50e-2 0.0133

Pretty up the table above for publication

conciseImmuneAlphaTable <- immuneAlphaTable %>% 
  dplyr::select(Type = type, Antigen = antigen, Month = month, Coef=Estimates, R2 = r2, P=`p.values`, Q = qval)
toAlphaTable <- conciseImmuneAlphaTable %>%
  mutate(
    Coef = cell_spec(format(Coef, digits = 3), "html",
                     bold = ifelse(P < 0.05, T, F),
                     italic = ifelse(P < 0.05 & Coef < 0, T, F),
                     background = ifelse(P < 0.05, ifelse(Coef < 0, "lightsalmon", "lightblue"), "")),
    R2 = cell_spec(round(R2, digits = 3), "html"),
    P = cell_spec(format_round(P, 3), "html",
                         bold = ifelse(P < 0.05, T, F),
                         background = ifelse(P < 0.05, 'yellow', '')
    ),
    Q = cell_spec(format_round(Q, 3), "html",
                         bold = ifelse(Q < 0.2, T, F),
                         background = ifelse(Q < 0.2, 'lightyellow', '')
    )
  ) %>%
   mutate(Antigen = gsub('ANY.ENV.PTEG', 'Any ENV PTEG', Antigen)) %>%
  mutate(Antigen = gsub('gp70_B.CaseA_V1_V2', 'gp70 B.CaseA V1-V2', Antigen))

toAlphaTable %>% kable("html", escape = F, digits = 3, align = 'c') %>%
  kable_styling("striped", "hover", full_width = F)  %>%
  collapse_rows(columns = 1:2, latex_hline = "full") %>%
  pass-> conciseAlphaTable.html

conciseAlphaTable.html %>% cat(file = 'tables/conciseAlphaTable.html')

conciseAlphaTable.html

immuneAlphaCompiledGaus <- use.immune %>% 
  group_by(type, antigen, month) %>%
  do(immuneValpha(., transformation = jac_box_cox_2)) %>% 
  ungroup

immuneAlphaTableGaus <- immuneAlphaCompiledGaus %>% mutate(qval = p.adjust(`p.values`, method = "BH"))
immuneAlphaTableGaus

## # A tibble: 20 x 8
##    type  antigen    month Estimates Standard.Errors p.values      r2  qval
##    <chr> <chr>      <dbl>     <dbl>           <dbl>    <dbl>   <dbl> <dbl>
##  1 CD4+  ANY.ENV.P…   6.5     30.5            12.6     0.016 1.34e-1 0.064
##  2 CD4+  ANY.ENV.P…  12       27.5            16.4     0.094 9.89e-4 0.235
##  3 IgA   gp41         0       -6.51           18.0     0.718 3.43e-2 0.818
##  4 IgA   gp41         6.5      2.15           15.5     0.89  6.04e-2 0.89 
##  5 IgA   gp41        12       32.9            16.4     0.045 8.62e-2 0.15 
##  6 IgA   p24          0       -6.77           18.5     0.714 2.82e-4 0.818
##  7 IgA   p24          6.5     15.9            17.9     0.374 3.31e-2 0.68 
##  8 IgA   p24         12       25.2            20.8     0.226 6.30e-2 0.452
##  9 IgG   Con.6.gp1…   6.5      6.48           17.2     0.706 2.44e-6 0.818
## 10 IgG   Con.6.gp1…  12        5.36           15.9     0.736 7.92e-5 0.818
## 11 IgG   gp41         0       10.1            15.4     0.51  1.19e-4 0.774
## 12 IgG   gp41         6.5     29.7             9.75    0.002 1.76e-1 0.01 
## 13 IgG   gp41        12       31.7             4.85    0     1.66e-1 0    
## 14 IgG   gp70_B.Ca…   6.5      3.24           18.5     0.861 2.06e-3 0.89 
## 15 IgG   gp70_B.Ca…  12       36.1             6.73    0     3.27e-2 0    
## 16 IgG   p24          0      -14.4            17.9     0.421 1.01e-1 0.702
## 17 IgG   p24          6.5     21.7            16.0     0.175 1.44e-1 0.389
## 18 IgG   p24         12       32.0            17.0     0.06  2.08e-1 0.171
## 19 IgG   ZM96.gp140   6.5      8.40           13.8     0.542 4.73e-4 0.774
## 20 IgG   ZM96.gp140  12       32.2             4.71    0     8.80e-2 0

conciseImmuneAlphaTableGaus <- immuneAlphaTableGaus %>% 
  dplyr::select(Type = type, Antigen = antigen, Month = month, Coef=Estimates, R2 = r2, P=`p.values`, Q = qval)
toAlphaTableGaus <- conciseImmuneAlphaTableGaus %>%
  mutate(
    Coef = cell_spec(format(Coef, digits = 3), "html",
                     bold = ifelse(P < 0.05, T, F),
                     italic = ifelse(P < 0.05 & Coef < 0, T, F),
                     background = ifelse(P < 0.05, ifelse(Coef < 0, "lightsalmon", "lightblue"), "")),
    R2 = cell_spec(ifelse(R2 < 0.01, format(R2, digits = 2),round(R2, digits = 2)) , "html"),
    P = cell_spec(format_round(P, 3), "html",
                         bold = ifelse(P < 0.05, T, F),
                         background = ifelse(P < 0.05, 'yellow', '')
    ),
    Q = cell_spec(format_round(Q, 3), "html",
                         bold = ifelse(Q < 0.2, T, F),
                         background = ifelse(Q < 0.2, 'lightyellow', '')
    )
  ) %>%
   mutate(Antigen = gsub('ANY.ENV.PTEG', 'Any ENV PTEG', Antigen)) %>%
  mutate(Antigen = gsub('gp70_B.CaseA_V1_V2', 'gp70 B.CaseA V1-V2', Antigen))

toAlphaTableGaus %>% kable("html", escape = F, digits = 3, align = 'c') %>%
  kable_styling("striped", "hover", full_width = F)  %>%
  collapse_rows(columns = 1:2, latex_hline = "full") %>%
  pass-> conciseAlphaTableGaus.html

conciseAlphaTableGaus.html %>% cat(file = 'tables/conciseAlphaTableGaus.html')

conciseAlphaTableGaus.html

psN1richP_gp41m6 <- psN1rich %>% sample_data %>% as.data.frame %>%
  ggplot(aes(x = MDS1, y = richEst, ymin = richLb, ymax = richUb, fill = as.factor(medcode_hl(IgG_gp41_Month_6.5)))) + geom_point(shape = 21, size = 4) + geom_errorbar() + scale_y_log10() + scale_fill_viridis_d(direction = -1, name = 'gp41 Primary Timepoint') + labs(x = "MDS1", y = "Richness (# ASVs)") + cowplot::theme_cowplot()
psN1richP_gp41m6

par <- options()
#options(repr.plot.width=6, repr.plot.height= 11)
#g <- grid.arrange(wuford_gp41, wuford_gp120, ncol = 2)
g <- cowplot::plot_grid(psN1richMDSP, psN1richP_gp41m6, ncol = 1, labels = c("A", "B"), label_size = 24, align = "v")
g

cowplot::save_plot('figures/richnessMDS1Gp41.png', g, base_width = 6, base_height = 8)

Next question. Do the donor and non donor groups ever have different immune responses?

This is in response to a reviewr comment. Strategy, write a function to, for one variable of interest, compare the groups, as I did for unifrac distance. I’ll use a simple logistic regression here. Similar to CapVar. I want a use.immune data set that includes whether people are a donor as a column, but that

immune.data %>%
mutate(isDonor = pub_id %in% muDoners) %>%
filter(
    (type == 'IgG' & 
    antigen %in% ants1 &
    month %in% c(6.5,12)
    ) |
    (type %in% c('IgG', 'IgA') &
     antigen %in% ants2 &
     month %in% c(0,6.5,12)
    ) |
    type == 'CD4+' &
    antigen == 'ANY.ENV.PTEG' &
    month %in% c(6.5, 12)
      )-> allparticipants.immune
head(allparticipants.immune)

## # A tibble: 6 x 14
##   visitno rx_code type  antigen   mag mag_bl response   day month ct   
##     <int> <chr>   <chr> <chr>   <dbl>  <dbl>    <int> <int> <dbl> <chr>
## 1       9 CTRL    IgA   gp41     177.   214.        0   182   6.5 CTRL 
## 2      12 CTRL    IgA   gp41     173    214.        0   364  12   CTRL 
## 3       9 CTRL    IgA   p24      966    454.        0   182   6.5 CTRL 
## 4      12 CTRL    IgA   p24     1419    454.        0   364  12   CTRL 
## 5       9 CTRL    IgG   Con.6.…    1      1         0   182   6.5 CTRL 
## 6      12 CTRL    IgG   Con.6.…    1      1         0   364  12   CTRL 
## # ... with 4 more variables: response_j <dbl>, assay <chr>, pub_id <dbl>,
## #   isDonor <lgl>

allparticipants.test <- allparticipants.immune %>% filter(visitno ==9, antigen == "Con.6.gp120.B",type == "IgG")
head(allparticipants.test)

## # A tibble: 6 x 14
##   visitno rx_code type  antigen    mag mag_bl response   day month ct   
##     <int> <chr>   <chr> <chr>    <dbl>  <dbl>    <int> <int> <dbl> <chr>
## 1       9 CTRL    IgG   Con.6.…     1    1           0   182   6.5 CTRL 
## 2       9 T2      IgG   Con.6.… 14938.   1           1   182   6.5 T    
## 3       9 T1      IgG   Con.6.… 30458.   1           1   182   6.5 T    
## 4       9 T4      IgG   Con.6.… 29971.   1           1   182   6.5 T    
## 5       9 T3      IgG   Con.6.… 31292.   7.75        1   182   6.5 T    
## 6       9 T4      IgG   Con.6.… 17209  362.          1   182   6.5 T    
## # ... with 4 more variables: response_j <dbl>, assay <chr>, pub_id <dbl>,
## #   isDonor <lgl>

donorTest <- function(x, transformation = medcode2, family = 'binomial'){
  loc_glm <- glm(transformation(mag) ~  isDonor, data = x, family = family)
  loc_glm %>% broom::tidy() %>% filter(term == 'isDonorTRUE') %>% dplyr::select(-term)
}
donorTest(allparticipants.test)

## # A tibble: 1 x 4
##   estimate std.error statistic p.value
##      <dbl>     <dbl>     <dbl>   <dbl>
## 1    0.380     0.506     0.751   0.453

Do on everything

allparticipants.immune %>%
  filter(ct == 'T') %>%
  group_by(type, antigen, visitno) %>%
  do(data.frame(donorTest(.))) %>%
  pass-> DonorEffectTable

DonorEffectTable

## # A tibble: 20 x 7
## # Groups:   type, antigen, visitno [20]
##    type  antigen            visitno  estimate std.error statistic p.value
##    <chr> <chr>                <int>     <dbl>     <dbl>     <dbl>   <dbl>
##  1 CD4+  ANY.ENV.PTEG             9  9.53e- 2     0.520  1.83e- 1   0.855
##  2 CD4+  ANY.ENV.PTEG            12 -1.17e-16     0.526 -2.22e-16   1.000
##  3 IgA   gp41                     2 -3.89e- 1     0.512 -7.60e- 1   0.447
##  4 IgA   gp41                     9 -1.72e- 1     0.518 -3.33e- 1   0.739
##  5 IgA   gp41                    12 -2.75e- 1     0.526 -5.23e- 1   0.601
##  6 IgA   p24                      2  1.29e- 1     0.509  2.54e- 1   0.799
##  7 IgA   p24                      9  3.65e- 1     0.521  7.00e- 1   0.484
##  8 IgA   p24                     12 -2.32e-16     0.524 -4.44e-16   1.000
##  9 IgG   Con.6.gp120.B            9  4.05e- 1     0.523  7.76e- 1   0.438
## 10 IgG   Con.6.gp120.B           12  1.33e- 1     0.516  2.58e- 1   0.797
## 11 IgG   gp41                     2  3.57e- 1     0.513  6.95e- 1   0.487
## 12 IgG   gp41                     9 -1.35e- 1     0.519 -2.59e- 1   0.795
## 13 IgG   gp41                    12  1.33e- 1     0.516  2.58e- 1   0.797
## 14 IgG   gp70_B.CaseA_V1_V2       9  4.05e- 1     0.523  7.76e- 1   0.438
## 15 IgG   gp70_B.CaseA_V1_V2      12 -1.33e- 1     0.516 -2.58e- 1   0.797
## 16 IgG   p24                      2 -1.64e- 1     0.510 -3.22e- 1   0.747
## 17 IgG   p24                      9 -3.92e- 2     0.528 -7.43e- 2   0.941
## 18 IgG   p24                     12 -4.01e- 1     0.520 -7.72e- 1   0.440
## 19 IgG   ZM96.gp140               9  9.53e- 2     0.520  1.83e- 1   0.855
## 20 IgG   ZM96.gp140              12 -4.42e- 1     0.521 -8.47e- 1   0.397

Ok. There appears to be no significant difference in magnitude group for any category.

allparticipants.immune %>%
  filter(ct == 'T') %>%
  group_by(type, antigen, visitno) %>%
  do(data.frame(donorTest(.,  transformation = function(x){jac_box_cox_2(x)},
                     family = 'gaussian'))) %>%
  pass-> DonorEffectGaus

Monkeys, I guess its debug aclock.

DonorEffectGaus

## # A tibble: 20 x 7
## # Groups:   type, antigen, visitno [20]
##    type  antigen            visitno estimate std.error statistic p.value
##    <chr> <chr>                <int>    <dbl>     <dbl>     <dbl>   <dbl>
##  1 CD4+  ANY.ENV.PTEG             9  0.441       0.256    1.72    0.0897
##  2 CD4+  ANY.ENV.PTEG            12  0.0673      0.265    0.254   0.800 
##  3 IgA   gp41                     2 -0.105       0.255   -0.413   0.681 
##  4 IgA   gp41                     9  0.236       0.259    0.911   0.366 
##  5 IgA   gp41                    12  0.349       0.260    1.34    0.184 
##  6 IgA   p24                      2  0.392       0.252    1.56    0.124 
##  7 IgA   p24                      9  0.274       0.258    1.06    0.293 
##  8 IgA   p24                     12  0.313       0.261    1.20    0.234 
##  9 IgG   Con.6.gp120.B            9 -0.0707      0.261   -0.271   0.787 
## 10 IgG   Con.6.gp120.B           12  0.0298      0.260    0.115   0.909 
## 11 IgG   gp41                     2  0.268       0.254    1.05    0.295 
## 12 IgG   gp41                     9 -0.113       0.261   -0.435   0.665 
## 13 IgG   gp41                    12 -0.0472      0.260   -0.182   0.856 
## 14 IgG   gp70_B.CaseA_V1_V2       9  0.185       0.260    0.710   0.480 
## 15 IgG   gp70_B.CaseA_V1_V2      12 -0.207       0.258   -0.799   0.427 
## 16 IgG   p24                      2  0.114       0.256    0.447   0.656 
## 17 IgG   p24                      9  0.242       0.264    0.916   0.363 
## 18 IgG   p24                     12  0.00988     0.260    0.0380  0.970 
## 19 IgG   ZM96.gp140               9  0.0605      0.262    0.231   0.818 
## 20 IgG   ZM96.gp140              12 -0.115       0.260   -0.442   0.660

Same deal when I do a normal logistic regression.

#save.image(file = "workspace.Rdata")